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Abstract. While a featureless, nearly scale invariant, primordial scalar power 
spectrum fits the most recent Cosmic Microwave Background (CMB) data rather well, 
certain features in the spectrum are known to lead to a better fit to the data (although, 
the statistical significance of such results remains an open issue) . In the inflationary 
scenario, one or more periods of deviations from slow roll are necessary in order to 
generate features in the scalar perturbation spectrum. Over the last couple of years, 
it has been recognized that such deviations from slow roll inflation can also result 
in reasonably large non-Gaussianities. The Starobinsky model involves the canonical 
scalar field and consists of a linear inflaton potential with a sudden change in the slope. 
The change in the slope causes a brief period of departure from slow roll which, in turn, 
results in a sharp rise in power, along with a burst of oscillations in the scalar spectrum 
for modes that leave the Hubble radius just before and during the period of fast roll. 
The hallmark of the Starobinsky model is that it allows the scalar power spectrum to 
be evaluated analytically in terms of the three parameters that describe the model, 
viz. the two slopes that describe the potential on either side of the discontinuity and 
the Hubble scale at the time when the field crosses the discontinuity. In this work, we 
evaluate the bi-spectrum of the scalar perturbations in the Starobinsky model in the 
equilateral limit. Remarkably, we find that, just as the power spectrum, all the different 
contributions to the the bi-spectrum too can be evaluated completely analytically and 
expressed in terms of the three paramaters that describe the model. We show that 
the quantity / NL , which characterizes the extent of non-Gaussianity, can be expressed 
purely in terms of the ratio of the two slopes on either side of the discontinuity in 
the potential. Further, we find that, for certain values of the parameters, / NL in the 
Starobinsky model can be as large as the mean value that has been arrived at from 
the analysis of the recent CMB data. We also demonstrate that the usual hierarchy of 
contributions to the bi-spectrum can be altered for certain values of the parameters. 
Altogether, we find that the Starobinsky model represents a unique scenario wherein, 
even when the slow roll conditions are violated, the background, the perturbations as 
well as the corresponding two and three point correlation functions can be evaluated 
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completely analytically. As a consequence, the Starobinsky model can also be used to 
calibrate numerical codes aimed at computing the non-Gaussianities. 



PACS numbers: 98.80.Cq, 98.70.Vc, 98.80.Es 



The scalar bi-spectrum in the Starobinsky model: The equilateral case 



3 



1. Introduction 

A nearly scale invariant primordial scalar power spectrum, along with the assumption of 
the concordant, background cosmological model (i.e. the spatially flat, ACDM model), 
seems to be in remarkable agreement with the data from the Wilkinson Microwave 
Anisotropy Probe (WMAP) of the anisotropies in the Cosmic Microwave Background 
(CMB) (for the most recent observations, see Refs. [1]; for earlier results, see Refs. [2]; 
for constraints on different inflationary scenarios, see Refs. [3, 4, 5, 6]). But, specific 
features in the primordial spectrum seem to improve the fit to the data to a reasonable 
extent (for an inherently incomplete list, see Refs. [7, 8, 9, 10, 11, 12, 13]), although the 
statistical relevance of such results remains difficult to assess. Notably, while a sharp 
drop in power on scales corresponding to the Hubble scale today has been found to fit the 
low CMB quadrupole better [7, 8, 9, 10], a burst of oscillations of a suitable amplitude 
and over a certain range of scales seem to provide a considerably improved fit to the 
outliers in the CMB angular power spectrum near the multipoles of I — 22 and 40 [11]. 
Moreover, interestingly, oscillating inflaton potentials that generate small modulations 
over a wide range of scales in the perturbation spectrum have also been found to perform 
better against the data than the more conventional power law spectrum [13]. 

As is well known, a featureless, scale invariant perturbation spectrum can be 
produced by a suitably long epoch of slow roll inflation (see any of the following texts [14] 
or reviews [15]). However, generating features in the scalar power spectrum require one 
or more periods of deviation from slow roll inflation [16, 17, 18]. For instance, it is found 
that one needs a large deviation from slow roll in order to produce a sharp drop in power, 
say, so as to fit the low quadrupole [16]. The larger the deviation from slow roll, the 
sharper the drop is found to be and, in fact, a brief departure from inflation — which 
leads to a sharp rise in power and a couple of oscillations before the spectrum turns 
nearly scale invariant — has been found to provide a good fit to the data at the lower 
multipoles [10]. In contrast, a small and short burst of deviation from slow roll seems 
sufficient to generate the oscillations that result in a much better fit to the outliers near 
the multipoles of I = 22 and 40 [11]. While certain features indeed lead to a better fit to 
the data than the standard power law primordial spectrum, the improvement in the fit 
is often arrived at with the introduction of extra parameters. As we pointed out above, 
the statistical significance of features achieved at the cost of additional parameters, say, 
from the Bayesian point of view [19], remains to be investigated satisfactorily. 

During the last few years, a variety of approaches have been developed to determine 
the extent of non-Gaussianities in the WMAP data (for an inexhaustive list, see 
Refs. [20]; in this context, also see the recent reviews [21]). Many of these analyses 
seem to indicate that the CMB may possibly possess reasonably large amount of non- 
Gaussianities. For instance, the WMAP seven year data constrains the parameter 
/ NL that is often introduced to characterize the extent of the non-Gaussianity to be 
/ NL = 32 ± 21 in the local limit, at 68% confidence level (see Ref. [1]; also see 
the reviews [21]). While a Gaussian primordial spectrum (which corresponds to a 
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vanishing / NL ) lies within 2-cr, evidently, the mean value seems to indicate a rather 
large amount of non-Gaussianity. 

The ongoing Planck mission [22] is expected to reduce the above-mentioned 
uncertainty in the WMAP's determination of / NL by a factor of about four or so. On 
the theoretical side, it has been realized that, if Planck indeed detects a reasonable 
amount of non-Gaussianity, then it can act as a powerful tool to substantially constrain 
the plethora of inflationary models that seem to be allowed by the currently available 
data [21]. For example, it has been established that the canonical scalar field models 
which lead to a nearly scale invariant primordial spectrum contain only a negligible 
amount of non-Gaussianity, with the magnitude of / NL turning out to be much smaller 
than unity [23, 24]. Hence, these models will cease to be viable if non-Gaussianity 
turns out to be substantial. However, it is known that primordial spectra with features, 
generated due to one or more deviations from slow roll inflation, can lead to reasonably 
large non-Gaussianities [25]. (It should be mentioned that, initial states other than the 
Bunch-Davies vacuum can also lead to features and, possibly, relatively substantial non- 
Gaussianities [26]. But, we shall not consider such situations in this paper.) Therefore, 
if non-Gaussianity indeed proves to be large, then, either one has to reconcile with the 
fact that the primordial spectrum contains features or one has to pay more attention to 
non-canonical scalar field models such as, say, the D-brane inflation models [27, 28]. 

Often, the efforts in evaluating the non-Gaussianities when features arise in 
the primordial spectrum due to a departure from slow roll have involved numerical 
computations [25], and it is instructive to consider a model wherein it may be possible 
to evaluate them analytically. In this work, we shall consider the Starobinsky model [16], 
which involves the canonical scalar field and consists of a linear inflaton potential with 
an abrupt change in the slope at a given point. The change in the slope causes a brief 
period of fast roll, which leads to a sharp rise in power and oscillations in the scalar 
spectrum for modes that leave the Hubble radius just before and during the period 
of fast roll. The important aspect of the Starobinsky model is that the scalar power 
spectrum can be evaluated analytically and, under certain conditions on the parameters, 
the analytic spectrum matches the exact spectrum, computed numerically, extremely 
well. We shall focus on the equilateral limit and illustrate that, remarkably, under the 
same conditions, all the contributions to the scalar bi-spectrum too can be calculated 
completely analytically and expressed entirely in terms of the three parameters that 
describe the Starobinsky model. We shall also discuss the ranges of the parameters over 
which the non-Gaussianity in the model can be large. 

At this point, we ought to mention that non-Gaussianities in the Starobinsky model 
have, in fact, been computed earlier [29]. But, since the formalism followed in the earlier 
work is rather different from the approach that we shall adopt, carrying out an exact 
comparison between the two efforts turns out to be difficult. While the non-Gaussianity 
parameter in the equilateral limit obtained in both the approaches lead to a sharp 
peak near the characteristic scale associated with the problem, we find the approach 
we have adopted here is able to capture the finer details as well. We believed that this 
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can attributed to the 5N formalism adopted in the earlier work, which is essentially 
applicable on super Hubble scales. Despite the additional efforts that were made to 
take into account the contributions due to the decaying mode — which become important 
when deviations from slow roll occur (in this context, see, for example, Refs. [30]) — we 
find that the approach still proves to be insufficient to arrive at the exact form for the 
final result. 

This paper is organized as follows. In the following section, we shall outline the 
essential aspects of the Starobinsky model. We shall first describe the background 
evolution in the model, and then go on to discuss the scalar power spectrum that arises 
in the model. In particular, we shall highlight the assumptions and approximations 
that are made in arriving at analytic expressions for the background evolution and 
the scalar power spectrum. In Sec. 3, we shall rapidly sketch the by-now standard 
procedure for evaluating the inflationary bi-spectrum of the scalar perturbations and 
the non-Gaussianity parameter / NL . In Sec. 4, using the method outlined in Sec. 3, we 
shall evaluate the dominant contribution to the scalar bi-spectrum in the Starobinsky 
model in the equilateral limit. In Sec. 5, we shall evaluate all the other sub-dominant 
contributions as well. In Sec. 6, we shall discuss the range of values of the parameters 
over which the non-Gaussianity parameter / NL can as large as indicated by the currently 
observed mean values. We shall also touch upon an issue related to the hierarchy of 
the different contributions to the bi-spectrum. Finally, we shall close in Sec. 7 with a 
brief summary and outlook. We shall relegate some details concerning the evaluation of 
certain integrals to an appendix. 

A few words on the conventions and notations that we shall adopt are in order 
at this stage of our discussion. We shall work in units such that h — c — 1, and 
we shall set the Planck mass to be M P1 = (8-kG) 1 ^ 2 . Also, while Greek indices shall 
refer to the spacetime coordinates, the Latin indices (barring the sub-script k which 
shall represent the wavenumber of the perturbations) shall denote the three spatial 
coordinates. Moreover, we shall work with the metric signature of (—,+,+,-1-). We 
shall express the various quantities in terms of either the cosmic time t, the conformal 
time f] (also, at times, denoted as r), or the number of e-folds N, as is convenient. As 
usual, while an overdot shall denote differentiation with respect the cosmic time, an 
overprime shall denote differentiation with respect to the conformal time coordinate. 
Lastly, a plus or a minus sign in the sub-script or the super-script of any quantity shall 
denote its value before and after the field crosses the discontinuity in the potential, 
respectively. 

2. The Starobinsky model 

In this section, we shall sketch the Starobinsky model [16] in some detail. We shall 
first describe the assumptions and the approximations that go into evaluating the 
background. We shall then go on to discuss the scalar modes and the resulting power 
spectrum that arise in the model. 
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2.1. The background evolution in the Starobinsky model 



Consider a spatially flat Friedmann universe that is described by the scale factor a, 
which is driven by the canonical scalar field 0. If the scalar field is described by 
the potential V(<j)), then the evolution of such a system is governed by the following 
Friedman and Klein-Gordon equations: 



H 2 



1 



3M 



+ V(<i>) 



and 



+ ZH<I> + Va = 0, 



(1) 



v(<f>) 



where H = a/a is the Hubble parameter, and = dV/d(p. As we had mentioned, the 
Starobinsky model consists of a linear potential with a sharp change in its slope at a 
given point. Let the value of the scalar field where the slope changes abruptly be 0o, 
and let the slope of the potential above and below O be A + and A_ , respectively. Also, 
let the value of the potential at = 0o be, say, Vq. In other words, the potential V(4>) 
describing the canonical scalar field is given by 

V + A + (0 - 0o ) for > 0o, , . 

V + A_ (0 - O ) for < O . 

This potential is displayed in the top left panel of Fig. 1. It is useful to note here that 
the quantities A + and A_ are of dimension three. 

Now, consider the case wherein the field is rolling down the above potential from 
a value such that > O . The exact trajectory of the field, computed numerically, is 
illustrated in Fig. 1 (top right panel, solid black line). In the slow roll approximation, 
the number of e-folds, TV = In (a/a{), is given by 

1 v 



N 



Ml 



14 



(3) 



where 0; and a ; are the initial values (i.e. at N — 0) of the field and the scale factor, 
respectively. If the field begins to roll down slowly, the trajectory of the field until it 
reaches the discontinuity at 0o can be easily obtained by carrying out the above integral 
for N. One obtains that 



0+ 



00 



0i - 00 + 



Vb 
A+ 



2 Mpj N 



1/2 



(4) 



and, as can been seen in Fig. 1 (top right panel, dotted blue line), this expression turns 
out to be an excellent approximation to the actual numerical solution (as one would 
expect, since the corresponding slow roll parameters are very small; see discussion 
below). This expression allows us to estimate the e-fold, say, iVo, at which the field 
reaches the point where the slope of the potential changes, and it is given by 



Nn = 



0i ~ 0C 

2MI 



00 + 



2Vo 

A + 



(5) 



The 'velocity' of the field before it reaches the discontinuity can be evaluated to be 

r 2 i - 1 / 2 

d0_u / , Vq 



diV 



-Mi 



0i - 00 + 



A, 



- 2 M p 2 , N 



(6) 
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Figure 1. The potential, the evolution of the scalar field and the first two slow 
roll parameters in the Starobinsky model. Top left panel: The potential (2) with 
4>o/M P1 ~ 0.707, Vq/M* = M i /M^ l ~ 2.37 x 10" 12 , A + /M|, ~ 3.35 x 10~ 14 and 
A-/Ml y ~ 7.26 x 10~ 15 , which corresponds to AA/A- = (A- - A + )/A- ~ -3.61. 
The mass M has been chosen such that the model is COBE normalized. Top right 
panel: The exact, numerically computed, evolution of the field cf> (for the above- 
mentioned values of the parameters) has been plotted as a function of the number 
of e-folds (the solid black line). The initial value of the field has been chosen to be 
</>i/M P1 ~ 0.849, and it corresponds to Nq — 10, where iVo is the e-fold at which the 
field crosses the discontinuity [cf. Eq. (5)]. The dotted blue line corresponds to the 
slow roll expression (4) for the evolution of the field. Clearly, the slow roll result is 
an excellent approximation to the actual result before the discontinuity. The dotted 
red line corresponds to Eq. (9), and it too accurately matches the exact solution after 
the discontinuity. Bottom left panel: The numerical result for the evolution of the 
first slow roll parameter t\ (the solid blue curve). The parameters that we work with 
lead to ei + ~ 10~ 4 [cf. Eq. (13)], as is confirmed by the plot. The dotted red curve 
represents the approximate expression of e\ after the discontinuity, viz. Eq. (14). 
Bottom right panel: The solid blue curve represents the exact evolution of the second 
slow roll parameter e2, obtained numerically. The dotted red curve represents the 
expression (17), which is valid after the field has crossed the discontinuity in the 
potential. The inset highlights the excellent agreement between the numerical and 
the analytical results. 
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Actually, the Starobinsky model assumes that the constant Vq is the dominant term 
in the potential for a range of <fi near O - This regime is very similar to the so-called 
vacuum dominated regime in hybrid inflation. In such a case, we can set V ~ Vq in the 
slow roll trajectory, and the expressions for <p + and d0 + /diV above simplify to 

A + Mp, AT , do. .1 . .1/-', 

The fact that, before the discontinuity, the slow roll trajectory can be well approximated 
by a straight line with a negative slope is also evident from Fig. 1 (top right panel). 
Further, in this limit, the above relation for iVo reduces to iVo = [Vo (<f>i — <po)/A + MpJ. 

When the field crosses <po, the slow roll approximation ceases to be valid, and 
we must return to the exact Klein-Gordon equation to understand the evolution of the 
field. Upon treating the number of e-folds as the independent variable, the Klein-Gordon 
equation in (1) can be rewritten as 

where e\ = —H/H 2 is the first slow roll parameter. Since the constant term Vq in the 
potential is dominant, as a first step, it is a good approximation to assume that the 
Hubble parameter H is constant. In other words, one always has e\ <C 1, even when the 
field passes over the discontinuity in the slope of the potential. This, in turn, implies that 
inflation does not end even though the slow roll approximation breaks down temporarily 
with the higher order slow roll parameters becoming large (see discussion below). In 
fact, this situation can be viewed as the very definition of inflation with a transition. 
The above Klein-Gordon equation can be easily solved under these conditions and, upon 
requiring the field and its derivative to be continuous at the transition, we obtain the 
following solution for the field when < <p Q : 



1 _ e -3(N-N ) 



3Hl 



(N-N ), (9) 



where 

AA = A_ - A + (10) 

and we have set Hq ~ Vb/(3Mpj). We have plotted this solution for the evolution of 
the field in Fig. 1 (top right panel, dotted red line). As can be noticed in the plot, the 
agreement with the exact, numerical, solution is very good. 

Let us now evaluate the slow roll parameters on either side of the transition due to 
the discontinuity in the slope of the potential. We have already introduced the first slow 
roll parameter e\. As is well-known, the second and the higher slow roll parameters are 
defined as [31, 32] 

dln l e ™l r m\ 
£n+i = dN for n > I. (11) 
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When the slow roll approximation is valid, the first two slow roll parameters can be 
expressed in terms of the potential V((f>) and its derivatives as follows: 



Ml 



ei 



v 



and e 2 ~ 2 Ml 



v 



V 



;i2) 



where = d 2 V/d(j) 2 . Clearly, these relations can be used until the field reaches O and 
one obtains that 



ei+ - 



M 2 



A. 



V + A + {4>- O ) 



A 2 



18 M*H$' 



(13) 



where the last expression has been arrived at upon assuming that V is dominant. The 
fact that the first slow roll parameter is almost a constant before the transition can 
be easily seen in Fig. 1 (bottom left panel). Further, since the potential is linear in 0, 
we have e 2+ — 4e 1+ , i.e. the second slow roll parameter is also a constant until the 
time the field reaches the discontinuity at 0o- Again, this behavior is confirmed by the 
corresponding numerical result plotted in Fig. 1 (bottom right panel). 

Once the field has crossed the break in the potential, one has to use the exact 
expressions to arrive at the behavior of the slow roll parameters. Upon using the 
solution (9), we find that the first slow roll parameter can be obtained to be 



AA 



1_ 2M 2 1 \dNj 18 MlHlY A. J' K ) 

This expression proves to be an excellent approximation to the exact result, as can be 
checked in Fig. 1 (it is represented by the dotted red curve in the bottom left panel). 
Let us now turn to the evolution of the second slow roll parameter after the transition. 
It can be written as 

which, we should emphasize, is an exact expression. In fact, we can arrive at the form 
of the second slow roll parameter after the transition, upon using the expression (14) 
for ei_ in Eq. (11). We find that, e 2 _ is given by 
_ 6AA e -3(v-v o) 
C2 ~ ~ ~A~ 1 - (Ka/aZ) e- 3 ^-^)- (16 ^ 
However, at sufficiently late times after the transition, when the slow roll evolution has 
been restored, we expect that e 2 _ ~ 4ei_, as it was in the case before the transition 
to fast roll (since the potential continues to be linear). Careful analysis points to the 
fact that, in order to achieve this relation at late times, one needs to actually take into 
account the property that, strictly, the Hubble parameter is not a constant. In fact, 
this amounts to using the usual approximation that consists of replacing the potential 
V by M 4 while its derivative is still calculated exactly. Upon taking this fact into 
account, one actually obtains that 



6AA 



-3(N-N ) 



£2- — 



A_ 1 - (AA/A-) e 



-3 (N-N ) 



4 6!. 



(17) 
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Figure 2. The solid blue curve denotes the evolution of the quantity de2/d(Ht), 
computed numerically, and plotted as a function of the number of e-folds. Note that 
we have worked with the same set of parameters as in Fig. 1. The dotted red curve 
represents the analytical expression (19). As is clear from the inset, the analytical 
expression approximates the actual numerical result very well after the transition. 
The vertical blue line is just a numerical artifact of the parametrization used here and 
we have carefully checked that it plays no role in the following considerations. 



with ei_ in the second term being given by Eq. (14). Clearly, as the first term vanishes 
at late times (i.e. at large N), the above expression satisfies the required relation 
between the first two slow roll parameters when slow roll has been restored much after 
the transition. Yet again, we find that this relation almost exactly matches the exact 
numerical result for 62, as is evident from Fig. 1 (see the bottom right panel). In 
particular, the above expressions point to the fact that the value of the second slow roll 
parameter soon after the transition is controlled by only one quantity, viz. AA/A_, and 
is given by 6 (AA/A-)/(l — AA/A-). Therefore, if |AA/y4_| is large, then €2 quickly 
approaches —6 immediately after the transition, as can also be checked in Fig. 1 (the 
figure actually corresponds to \AA/A-\ — 3.61). 

As we shall discuss in the next section, in order to evaluate the dominant 
contribution to the scalar bi-spectrum in the Starobinsky model, we shall also require 
the quantity 62- One can show that, it can be written as 

^ = + l2He, -3He 2 -AHe\ + 5He 1 e 2 - |e^, (18) 
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and, under the conditions of our interest, this relation reduces to 

18i/ AA e -HN-N ) 

— — (19) 

A- [1 - (AA/A_) e - 3 ( N - N o)] 2 ' 

We should mention here that, in arriving at the final equality, we have ignored the 
second term involving ei_ in Eq. (17). The fact that this is a valid approximation is 
confirmed by the numerical analysis. As in the case of e 2 _, we find that, during the 
period soon after the transition, the value of e 2 _ is determined only by ratio AA/A-. 
In Fig. 2, we have plotted the above expression for e 2 as well as the exact numerical 
result. It is clear from the figure that the expression (19) is an excellent approximation 
to the actual result. 

Let us now briefly pause to summarize the results we have obtained. It is evident 
from the above expressions for the slow roll parameters that the change in the slope 
in the potential leads to a departure from slow roll for a short span of time, before 
slow roll is again recovered when the friction on the scalar field due to the expansion of 
the universe begins to dominate the force due to the potential. While the first slow roll 
parameter e\ remains small throughout the evolution (which essentially occurs due to V 
being the dominant term in the potential around 4>o), the second slow roll parameter e 2 
and its time derivative e 2 exhibit a sharp rise as the field crosses <p , before falling off. 
We have been able to obtain simple analytical expressions for the various background 
quantities even in the regime where the conditions for slow roll are violated. This 
point turns out to be the key for our analysis. As we shall see, the simple expressions 
for the slow roll parameters permit us to evaluate the bi-spectrum in the equilateral 
limit completely analytically. However, before we turn to the evaluation of the scalar 
bi-spectrum, let us discuss the power spectrum that arises in the Starobinsky model. 

2.2. The scalar power spectrum in the Starobinsky model 

Let TZ(r], x) be the dimensionless curvature perturbation induced by the scalar field, and 
let TZk(fj) denote the associated Fourier modes defined through the relation 

/d 3 k 
(2^F^ ke * k ' x (20) 

On quantization, the operator corresponding to the curvature perturbation can be 
expressed as 

= J j0j- 2 [a^h^e^ + atf^e-^] , (21) 

where and a t are the usual creation and annihilation operators which satisfy the 
following non-trivial commutation relation: [ak, aj,] = (k — p). In the case of the 
canonical scalar field, the modes f k are governed by the differential equation [14, 15] 

f> + 2-f k + k 2 f k = ^ (22) 

z 
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where z = a<p/ H = aM Pl y/2 e±. The quantity TZ^ is of dimension —3, fk of dimension 
—3/2, while z is of dimension one. It is useful to note that, in terms of the Mukhanov- 
Sasaki variable Vk = z fk (which is of dimension —1/2), the above equation for reduces 
to 



z" 



v'i + [k 2 --j Vk = 0. (23) 

The dimensionless scalar power spectrum V s (k) is defined in terms of the correlation 
function of the Fourier modes of the curvature perturbation as follows: 

(OIKMKMIO) = ^ V s (k) (k + p) , (24) 

where the vacuum state |0) is defined as ak|0) = 0, V k. Since we can write, 
= (Sk fk + a^k /it)' one obtains that 

".<*) = £ iai' = £ (¥Y (25) 



with the right hand side (which depends on time) evaluated in the super Hubble limit 
[i.e. when k/(aH) — > 0] or, more generically, at the end of inflation. The curvature 
perturbation is usually assumed to be in the Bunch-Davies vacuum, which corresponds 
to choosing Vk — > e lkri /\/2~k in the sub Hubble limit, i.e. as k/(aH) — > oo. It is clear 
from Eq. (23) that it is the 'effective potential' z" / z which determines the evolution 
of the scalar perturbations. It can be written in terms of the slow roll parameters as 
follows: 

z" „, 2 / 3 1 2 1 1 \ 

_ = n 2 i 2 - ei + - e 2 + - e\ - - Cl e 2 + - e 2 e 3 J , (26) 

where % = a' /a = a H is the conformal Hubble parameter, and it should be emphasized 
that the above expression is an exact one. 

Our aim now is to derive the scalar power spectrum for modes that leave the 
Hubble radius in the vicinity of the transition at 0o- Before the transition to a brief 
period of fast roll near <p Q , all the slow roll parameters remain small, and the effective 
potential simplifies to z"/z ~ 2H 2 . We should emphasize here that, for simplicity, 
we are neglecting terms involving the first order slow roll parameters that would lead 
to a non-vanishing spectral index for modes far from the characteristic scale, viz. the 
mode which leaves the Hubble radius when the scalar field crosses the break in the 
inflaton potential. Recall that, since the evolution is dominated by the constant term 
Vo in the potential, the expansion is of the de Sitter form corresponding to the Hubble 
parameter H . Therefore, around the transition, the scale factor can be expressed in 
terms of the conformal time as a(i]) = —(Hot])' 1 , so that we have z" j z ~ 2 / r/ 2 . In 
such a case, the solution to the Mukhanov- Sasaki variable Vk that satisfies the standard 
Bunch-Davies initial condition is given by 

= 4^(^-1^] e~ lkv - (27) 



2k V kl l 
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During the transitory fast roll regime, as the field crosses 0o, the slow roll parameters 
are no longer small and, hence the solutions to the Mukhanov-Sasaki variable v k can 
be expected to look different. In particular, we have seen that the second slow roll 
parameter €2 and its time derivative €2 can be large immediately after the transition. 
From Eq. (26), it is then clear that, a priori, one can no longer expect the effective 
potential to be just given by z"/z ~ 2'H 2 . However, remember that e\ remains small 
even during the transition. Further, one has e 2 = He 2 e 3 . Upon using these facts and 
the result that e 2 can be approximated using Eq. (18) (and, of course, the property that 
V<t>4> — 0), one finds that, after the transition 

Z" n.2 f n 3 1 2 1 

~z^ U 2+ 2 e2 - + 4 e - + 2 e2 - e3 - 



U 2 



3 1 o 1 / _ 1 



2+2 6 2 _ + i 6 2 _ + - l-3e 2 _--e 2 _ 



2 



2'H 2 . (28) 



This implies that certain cancellations occur so that the effective potential still retains 
the same form, i.e. z"/z ~ 2'H 2 , even after the field has crossed 0o- But, due to the 
fast roll, post-transition, the modes v k do not remain in the Bunch-Davies vacuum and, 
as a result, the solution to v k takes the general form 

7^(77) = -^=(1-^-) e"^ + J*= ( 1 + ±) e lfc ", (29) 
k V2k V kr lJ V2k \ krjj 

where a k and (3 k are the standard Bogoliubov coefficients. 

The Bogoliubov coefficients a k and f3 k can now be determined by carefully matching 

the mode v k and its derivative at the transition at <p . In order to do so, we first need 

to understand the behavior of z and the effective potential z" / z across the transition. 

In terms of the conformal time coordinate, the transition occurs at r] = — (a i/o) \ 

where ao is the scale factor at 770- Recalling that z = a M Pl y/2e[, and upon using the 

expression of e\ before and after the transition, we arrive at 

"I" alAA^ (30) 



3H3 V 

As a consequence, one obtains that 



for 11 < 770, 
for 77 > i] . 



S ffi rf for r, < r i0 , 



for 77 > 770, 



3# 3 77 2 3 
and it is clear that z' jumps at the transition by the amount 

n 2 A A 

W(m)h = z'Mo) -<(Ub) = — (32) 

-no 

So, at the transition, that is to say when 77 ~ 770 (and only at the transition!), we can 
write 

^^-^efo-,*), (33) 

o -no -no 
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where Q(x) is the step function. This implies that the effective potential z"/z can be 
described as a Dirac delta function at the transition, and we have 

z" 3 a H AA 



(77-770). 



(34) 



The above expression for z"/z allows us to establish the matching conditions on 
the modes at the transition, which read as 



[v k (Vo)}± = v k (r/o) - v+ (770) = 



and 



W k (m)}± = VM - wit' (770) = 



3 a H AA 



Vk (?7o) • 



(35) 



(36) 



The modes (27) and (29), along with these two conditions, then lead to the following 
expressions for the Bogoliubov coefficients a k and (3 k : 

3iAA k 



a k = 1 + 



2A+ k 
3iAA k 
2A+ ~k 



1 + 



k 2 
k 2 



1 + 



k 



3 2ife/feo 



(37) 
(38) 



where ko = — 1/770 = ao-ffo corresponds to the mode that leaves the Hubble radius at the 
transition. The ratio k/k = k/(a H ) can also be viewed as the ratio of the physical 
wavenumber at the transition, viz. k/ao, to the characteristic physical wave number 
^phys ^ Therefore, upon using the modes (29) post-transition, the resulting power 
spectrum, evaluated as k 77 — > 0, is found to be [16] 



r s (k) 



+ 



Ho 
2ir 

Ho 
2ir 

2ko 

k 
2ko 

k 



cos 



sin 



3^ 

A_ 

3HI 
A_ 
2k 
ko 
2k 
ko 



dk - /3fc| 2 

SAAko 
~ A + ~k 
9AA 2 k 2 
2A 2 k 2 



k 2 
k 2 



sin 



+ 



1 + 



kr 



1 - 



k 2 

k 2 



cos 



k 2 
2k 
k 



1 + 



2k 
k 

> 

6 



k? 



A- 



(39) 



and it should be emphasized that the spectrum depends on the wavenumber only 
through the ratio k/k . Also, as one would expect, this power spectrum turns scale 
invariant far away from k on either side. While, for k/ko — > 0, its scale invariant value 
is given by 



Jim n(*0 = 

k/k -> 

for k/ko — > 00 ; if simplifies to 
hm V s (k) = 

k/ko 00 



Ho 
2tt 



Ho 
2vr 



3HI 

A, 



3H 2 



A_ 



(40) 



(41) 
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Figure 3. The scalar power spectrum in the Starobinsky model. While the blue solid 
curve denotes the analytic spectrum (39), the red dots represent the corresponding 
numerical scalar power spectrum that has been obtained through an exact numerical 
integration of the background as well as the perturbations. The green dashed lines 
indicate the asymptotic scale invariant values [cf. Eqs. (40) and (41)]. And, the inset 
highlights the small difference between the exact numerical spectrum and the analytic 
one over a certain window in the wavenumber. In plotting these spectra, we have 
worked with the same values of the parameters that we had mentioned in the first 
figure. The first scale (i.e. the largest one) for which the power spectrum has been 
numerically computed corresponds to the physical wavenumber given by kja\ = 500 Hi 
at the beginning of inflation, which is an arbitrary but convenient choice. In terms of 
this scale, ko is given by k = yk{, where 500 Hi ye~ N " — H . Given that Nq = 10 and 
Hi ~ Hq (as can be easily checked in the slow roll approximation and as can be verified 
with great accuracy numerically), one obtains fc n ~ 44fc;. Evidently, the analytically 
evaluated spectrum is in remarkable agreement with the numerical result. We should 
mention that an inaccurate determination of ko would have resulted in a 'phase shift' 
between the numerical and analytical spectra. 



Given a H , or equivalently, Vq, COBE normalization (i.e. the amplitude of the power 
spectrum at suitably small scales) determines the value of A_. In Fig. 3, we have plotted 
the above analytical power spectrum and the corresponding numerical result for certain 
values of the parameters for which COBE normalization is achieved and the assumptions 
of the Starobinsky model are valid. The two results are clearly in good agreement, which 
indicates the remarkable extent of validity of the assumptions and approximations that 
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are used to arrive at analytic forms for the background evolution, the modes as well as 
the power spectrum in the Starobinsky model. 

3. A rapid outline of the procedure for evaluating the scalar bi-spectrum 

In this section, we shall quickly sketch the by-now commonly used procedure for 
evaluating the scalar bi-spectrum generated during inflation (for the original discussion, 
see Ref. [24]; for further discussions, see Refs. [25, 27] and also the recent review [28]). 
Due to the constrained nature of the Einstein equations, the scalar perturbations 
during inflation can be basically described by a single function, say, the curvature 
perturbation 1Z. Just as the scalar power spectrum can be expressed in terms of the 
two point function of the curvature perturbation [using Eqs. (21) and (24)], the scalar 
bi-spectrum is related to the three point function through a suitable Fourier transform 
[cf. Eqs. (68) and (70)]. It is well known that, at the linear order in the perturbations, 
the curvature perturbation TZ is governed by a quadratic action [15], which leads to 
the equation of motion (22) that we had considered in order arrive at the the power 
spectrum in the previous section. Evidently, in a linear and Gaussian theory, the three 
point function will be identically zero§. Therefore, in order to evaluate the bi-spectrum 
during inflation, the first step that needs to be taken is to arrive at the action describing 
the curvature perturbation at the next higher, i.e. the cubic, order. Then, based 
on the cubic order terms, one evaluates the corresponding three point function of the 
curvature perturbation (and, thence the bi-spectrum) using the standard techniques of 
perturbative quantum field theory. In what follows, we shall first describe as to how 
the cubic order action is arrived at, and then discuss the contributions to the scalar 
bi-spectrum due to the various terms in the cubic order action. 

3.1. The cubic order action describing the curvature perturbation 

The action that describes the curvature perturbation at the cubic order is usually arrived 
at using the Arnowitt-Deser-Misner (ADM) formalism [33]. Recall that, in the ADM 
formalism, the metric is expressed in terms of the lapse and the shift functions N and ./V* 
as follows: 

ds 2 = -N 2 d{x ) 2 + h i:j (N l dx° + dx*) (N 3 dx° + dx j ) , (42) 

where x° and x l denote the time and the spatial coordinates. As is common 
knowledge, the lapse and the shift functions N and N l turn out to be Lagrangian 
multipliers in the action for the complete system, and the variation of the action 
with respect to them leads to the Hamiltonian and the momentum constraints. The 
remaining equations govern the dynamics of the spatial metric hij. These six equations 
(corresponding to the independent components of the spatial metric) and the four 

§ It may be useful to point out here that, even a linear theory can lead to non-Gaussian distributions, 
when one works with, say, non-vacuum initial states. 
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constraint equations essentially constitute the ten Einstein equations corresponding to 
the complete spacetime metric. 

The system of our interest is gravity described by the Einstein-Hilbert action and 
a scalar field that is governed by the canonical action. For such a case, in terms of the 
metric (42), the action describing the complete system can be written as [24, 28] 



S[N, Nt, h i:j ,(l)} = J dx° J d 3 x N Vh 

d(j)\ 2 N l d(f) 



l -{E t] E^ -E 2 )+ ®R 



N 2 



+ 
+ 



2N 2 \dx° 
TV* N j 



2N 2 



N 2 dx° 
d4d j( t>-V{(t>) 



die 



-tiid i( f>d j( t> 



(43) 



where h = det (hij), is the spatial curvature associated with the metric hij, and 
the scalar field is assumed to be, in general, dependent on time as well as space. The 
quantity Eij is the rescaled second fundamental form and is given by 



E ij = NK ij = ± 



dhjj 
dx° 



((3) V . Nj + 0) V . N .j 



(3)i 



(44) 



while E = h^ E l i . Evidently, in the above action, the quantities within the first and 
the second square brackets correspond to the gravitational and the scalar field parts, 
respectively. 

In the Hamiltonian formulation that we are working in, the conjugate momenta 
corresponding to the spatial metric and the scalar field 0, are given by 



7T 



Ml Vh 



(K ij - ti j K) and vr 



Vh ( d(j) 



N \dx° 



(45) 



The Hamiltonian density, say, H, of the entire system can be expressed in terms of the 
generalized coordinates h^ and 0, and the corresponding conjugate momenta iiij and 
7r^ as follows: 



H = 




7T lJ 7T, 



71 

~2 



MlNVh 



Vh , ., 



+ _ h H di< p dj(f)+ VhV( ( j ) ))+ 7T N l d t( j), 



(46) 



where, clearly, the first line corresponds to gravity, while the second line corresponds 
to that of the scalar field. Moreover, it is evident from the above Hamiltonian that the 
lapse and the shift functions indeed act as Lagrange multipliers, as we mentioned above. 
Varying with respect to N and iVj respectively leads to the following Hamiltonian and 
momentum constraints: 



Ml 
2 

+ 



[(K lj - h ij K) (K^ - K)-2K 2 - (3) i?] 



2N 2 \dx° 



d<p 



N'dd) +-ti j d i <f>d j <j> + V( ( f>) = 0, 



(47) 
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and 

-2 (3) Vi7r ij + ir^ h ij 84 = 0. (48) 

Let us emphasize that, so far, the formalism has been general. Nothing has yet been 
said about the symmetries and the form of the metric tensor, nor has any approximation 
been made. 

In the case of the Friedmann metric, based on the time coordinate that one works 
with, the lapse function N can either be set to unity if one chooses x° to be the cosmic 
time t or we can set N = a when working with the conformal time x° = rj. At this stage, 
it is instructive to briefly compare the ADM approach with the more common method 
of arriving at the equations governing the perturbations from the Einstein equations. In 
the latter approach, when the scalar perturbations are taken into account, the spatially 
flat, Friedmann metric in an arbitrary gauge is usually written as 

ds 2 = - (1 + 2 if) dt 2 + 2 a(t) diB dt dx i + a 2 {t) \{l - 2 tjj) dx i dx j 

+ 2 didjS dx i dx j ] , (49) 

where tp, B, if) and £ are the scalar functions that describe the perturbations. Upon 
comparing this line-element with Eq. (42), it is clear that, we have 

Nirf-N 2 = -(l + 2<p), Ni = adiB, 

h,^ = a 2 (1 - 2 ip) Stj + 2 a 2 d^E. (50) 

At the first order in perturbations, these equations reduce to 

N ~ 1 + <p, N t ~a d t B, h^ ~ a 2 (% - 2 tjj 5 tj + 2 , (51) 

and the lapse function N is just unity at the zeroth order because we have chosen to 
work in terms of the cosmic time coordinate. 

From the above expressions for the different components of the metric tensor, one 
can evaluate the second fundamental force and the spatial curvature, and one obtains 
that 

Kij =(h~ ZHiP -ijj- H^j a 2 5ij - adidjB + a 2 didj (s + 2H£^j, (52) 

( 3) R = 4 di&rj). (53) 
a 2 

For a consistent perturbation theory, we must, of course, perturb the scalar field as 

well. Let us write the scalar field as a homogeneous term plus a perturbative part, 

i.e. as + 5(f). If we then use the above expressions for and in the constraint 
equations (47) and (48), we arrive at 

-3H (h<p + i£) + ^did i i)-—d i d i B + Hd i d i 8 
V /a 2 a 

= 2~W V + * ^ + V(t> ^) ' (54) 
0V + ff0V=2^f0Vty. (55) 
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These are nothing but the time-time and time-space perturbed Einstein equations. 

Our main goal now is to evaluate the action (43) at the cubic order in the perturbed 
quantities. Ideally, one would have liked to work in an arbitrary gauge, as is done to 
arrive at the standard action that describes the curvature perturbations at the quadratic 
order [15]. However, as one can imagine, the calculation in an arbitrary gauge turns 
out to be rather complicated and, to the best of our knowledge, so far, the action at 
the cubic order seems to have been arrived at by working in a specific gauge. In fact, 
the gauge that proves to be the most convenient is the one wherein £ = and 5<p = 0, 
i.e. the gauge wherein perturbations in the inflaton are assumed to vanish [24, 27, 28]. 
Then, using the perturbed equations derived above, it is easy to show that 

<P = -^, aB = ^- ° 2 f (56) 
H H 2MIH 2 K ' 

where the operator d~ 2 is defined by the relation d~ 2 did l tp = ip. The first of these two 
equations can be obtained from the momentum constraint, while the second arises due 
to the Hamiltonian constraint. Let us now introduce the curvature perturbation 71, a 
quantity that is conserved at super Hubble scales, and is given by [34, 35] 

7l = -iP-S-5<P, (57) 



In the gauge that we are working in, we have 71 = —ip, and the above equations simplify 
to 

These expressions match Eqs. (28) and (29) in the first article of Refs. [27], if one notices 
that, in the case of a canonical scalar field, the quantity E introduced in the article is 
given by £ = 2 /2. 

We shall assume that the spatial metric hij is given by 

h tj = a 2 (t)e 2n ^ S tj (59) 

and, it is worth noting that, in the gauge that we are working in, the quantity 7Z that 
appears in the exponential above is essentially the curvature perturbation. One finds 
that, it suffices to solve the two constraint equations for the lapse and the shift functions 
at the linear order in perturbations. This is due to the fact that, at the cubic order 
in 71, two type of terms arise in the action which involve the lapse and shift functions at 
the second and the third orders. However, it is found that, the coefficient of the second 
order term contains a first order constraint, while the third order term proves to be 
proportional to the zeroth order constraint, both of which, evidently, vanish [24, 27, 28]. 

For the canonical scalar field of our interest, the action that governs the curvature 
perturbation 7Z at the quadratic and the cubic orders can be arrived at from the original 
action (43). After a considerable amount of manipulations, most of which involve using 
the solutions to the constraint equations, repeatedly integrating the action by parts and 
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systematically throwing away the surface terms, one can show that, at the quadratic 
order in K, the action (43) simplifies to [24, 27, 28] 

1 



s 2 [n\ 



d?7 / d J x z 1 



n' 2 - {dnf 



(60) 



It is straightforward to check that the variation of this action with respect to 1Z leads 
to the differential equation (22) that we had considered earlier. 

The action at the cubic order in the curvature perturbation can be obtained in a 
similar fashion, and is found to be [24, 27, 28] 



S 3 [R] = M^J d v J d 3 



a 2 e 2 1 1Zn' 2 + a 2 e 2 i n(d1Zf 

2 



-2ae 1 K' (d'll) (d iX ) + y ci e' 2 1Z 2 K> + | {&K) (d iX ) {d\) 



+ 1 (d 2 n) {d x ) 2 



a J 7 



6C2 
51Z 



where the quantity X is defined through the relation 

X = d~ 2 A, 
with A being given by 



A 



1l = ae l 'R,'. 



(61) 



(62) 



(63) 



2 M 2 if 2 

Note that, while A is of dimension one, x is °f dimension —1. The quantity 
5C2/51Z denotes the variation of the Lagrangian density corresponding to the quadratic 
action (60), and can be written as 

5C 2 



511 



A + HA-e 1 (d 2 TZ). 



(64) 



The term J r (5£ 2 /^) that has been introduced in the above cubic order action refers 
to the following expression: 



5C2 
511 



1 



a 2 He 2 1Z 2 + 4a1Z1Z' + (d'ft) (d iX ) - — {Oil) 2 

H 



2aH 

+ [A{d t 1Z) + {d 2 1Z){d lX )\ 5 ij d : 



5C2 
51Z 

5C2 
51Z 



5_C2 
51Z 



(65) 



Before we proceed, four remarks are required to be made. Firstly, we should point 
out that, in the above expressions that constitute the action S3, the spatial indices 
are to be raised or lowered with the Kronecker symbol. Secondly, notice that, all the 
terms in the expression (65) for J r (5C2/51Z), barring the first one, contain a derivative 
of the curvature perturbation (either a time or a spatial derivatives or both). Since 
the quantities have to be evaluated at the end of inflation, the terms which involve the 
derivatives will not contribute to the final results. Thirdly, we have checked that our 
expression for S 3 is consistent with the results obtained earlier (such as in the first three 
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articles of Refs. [27]). Lastly, it is useful to note that, because the speed of sound is 
unity for the canonical scalar field, no term involving TZ 3 arises, as it happens in cases 
involving non-canonical scalar fields. 



3.2. The definition of the bi-spectrum and the different contributions 

The goal now is to treat the third order action (61) as the interaction term and evaluate 
the three point function of the curvature perturbation using the standard techniques of 
perturbative quantum field theory. To begin with, it can be shown that the last term 
in the action (61) which involves 5C2/S1Z can be removed by a field redefinition of 1Z of 
the following form: 

n^n n + F(7i n ), (66) 

where F = e 2 7?. 2 /4 (for further details, see Refs. [24, 27, 28]). With such a redefinition, 
the interaction Hamiltonian corresponding the action S 3 [1Z] above can be written in 
terms of the conformal time coordinate as 



H iQt ( V )= - Ml J d 3 x 



a 2 einn' 2 + a 2 t[n{dny - 2ae l n'(d i n)(d iX ) 



+ yde^tt' + |(^)(%) {d 2 X ) + j {d 2 n) (d X f 



(67) 



The three point correlation function of the curvature perturbation can be expressed 
in terms of the Fourier modes as follows: 

(^,x)^,x)^,x))= J J / ApL 

x (^(77)^(77)^3(77)) e^ kl+k2+k3 K (68) 

At the leading order in the perturbations, one then finds that the three point correlation 
in Fourier space is described by the integral [24, 27, 28] 

(7t kl (77 e )7t k2 (77 e )7t k3 (?7e)) 



Jili 



dr 



^k 1 (r/e)^ k2 (r/e)7t k3 (r/e),^int(r) ), (69) 



where H mt is the operator corresponding to the interaction Hamiltonian (67), while 771 
is the time at which the initial conditions are imposed on the modes when they are well 
inside the Hubble radius, and ?7e denotes a very late time, say, close to when inflation 
ends. Moreover, while the square brackets imply the commutation of the operators, 
the angular brackets denote the fact that the correlations are evaluated in the initial 
vacuum state (viz. the Bunch-Davies vacuum in the situation of our interest). 

The scalar bi-spectrum £> s (ki, k 2 , k 3 ) is related to the three point correlation 
function of the Fourier modes of the curvature perturbation, evaluated at the end of 
inflation, say, ?7 e , as follows [1]: 

(^ kl (r7e)^ k2 (r7e)^ k3 (r7e)) = (27r) 3 /3 s (k 1 ,k 2 ,k 3 )5( 3 )(k 1 + k 2 + k 3 ). (70) 
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For convenience, we shall set 

B^MM) = (2?rr 9/2 G(k u k 2 M) (71) 
and, it is useful to note that, the two quantities B s and G are of dimension —6. We can 
arrive at the quantity G(k 1 , k 2 , k 3 ) from the expressions (69) and (67) for the three point 
function in Fourier space and the interaction Hamiltonian Hi nt , respectively. Upon using 
the decomposition (21), as well as the Wick's theorem, which applies to the correlation 
functions of free quantum fields, we find that G(ki, k 2 , k 3 ) can be written as 

7 

c=i 



C=l 



+ [KM KM tiM] ^(ki,k 2 ,k 3 ) 

+ G 7 (k 1 ,k 2 ,k 3 ). ' (72) 

The quantities Q c (k 1 , k 2 , k 3 ) with C = (1,6), which are of dimension —7/2, correspond 
to the six terms in the interaction Hamiltonian (67), and are described by the 
integrals [24, 25, 27] 

G 1 (k 1 , k 2 , k 3 ) = 2 i / dr a 2 e\ (f* ki f* f» + two permutations) , (73) 
2 % (ki ■ k 2 + two permutations) 

dra 2 elr ki r k2 r k3 , (74) 
ki • k 2 



^ 2 (ki,k 2 ,k 3 ) 



x 



£ 3 (k 1; k 2 ,k 3 ) = -2i dra 2 ei 



+ five permutations 



/* ft* 
k\ J k2 J ^3 



(75) 



^ 4 (ki, k 2 , k 3 ) = i J dr a 2 e x e' 2 (f kl fl 2 f k \ + two permutations) , (76) 



g 5 (k u k 2 M) 



- f C dra 2 e\ 



ki • kr 



k 2 



fcl J ki J kz 



+ five permutations 



(77) 



dr a 2 e\ 



v k 2 (k 2 -k 3 

k 2 k 2 
/i 2 At 3 



ki J k2 J k$ 



+ two permutations > . 



(78) 



The additional, seventh term G 7 (ki, k 2 , k 3 ) arises due to the field redefinition (66), and 
its contribution to G(k 1 ,k 2: k 3 ) is found to be given by [24, 27, 28] 



G 7 (ki,k2,k3) 



{\fk 2 (Vc)\ 2 \fk 3 (Vc)\ 2 + two permutations) . (79) 
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Note that, whereas the first three integrals Qi, Q 2 and Q% involve e 2 , the next 
three, viz. Q 4 , Q 5 and Q e , depend on either e\ e' 2 or e\. And, evidently, the term G 7 is 
proportional to e 2 . During slow roll, the parameter t\ is almost constant and is typically 
of the order of 1CT 2 or so, while the second and higher slow roll corrections are even 
smaller than e\. Therefore, in a slow roll inflationary scenario driven by the canonical 
scalar field, it is found that the first three terms G\, G 2 and G 3 (that involve the 
integrals Qi, Q 2 and Q 3 ) and the term G 7 , which arises due to the field redefinition, that 
contribute the most to the bi-spectrum [24, 27, 28]. However, when there are deviations 
from slow roll, it has been noticed that the contribution due to the term G 4 proves to 
be the largest to the bi-spectrum as the corresponding integral Q± contains the quantity 
eie' 2 [25]. 

3.3. The non-Gaussianity parameter / NL 

The observationally relevant dimensionless non-Gaussianity parameter / NL that we had 
discussed about in the introductory section is related to the three point correlation 
function of the curvature perturbation as follows. It is introduced through the 
equation [24, 25] 

^,x)=^x)-^ [^,x)] 2 , (80) 

where 1Z G denotes the Gaussian quantity, and the factor of 3/5 arises due to the 
relation between the Bardeen potential and the curvature perturbation during the 
matter dominated epoch. In Fourier space, the above equation can be written as 

Kt = K c_^|_^ K c K c_ p (81) 

Upon using this relation, the Wick's theorem, which applies to the correlation functions 
of free quantum fields or, equivalently, to Gaussian random fluctuations, and the 
definition (24) of the power spectrum, one can arrive at the three point correlation 
of the curvature perturbation in Fourier space in terms of the parameter / NL . It is found 
to be 

x [kf V s (k 2 ) V s (k 3 ) + two permutations] . (82) 

Using this expression for the three point function and the definition (70) of the bi- 
spectrum, we can then arrive at the following relation between the non-Gaussianity 
parameter / NL and the bi-spectrum B s (ki, k 2 , k 3 ): 

1 

/ NL = - - (2 tt)- 4 (2 tt) 9 / 2 kf kl kl B s (kx, k 2 , k 3 ) 
x [^1 'Psikz) ^si^) + ^ wo permutations] 1 
= -^(2 7T)- 4 ^A; 2 3 A;3 3 G(k 1 ,k 2 ,k3) 

x [kf V s (k 2 ) V s (k 3 ) + two permutations] _i . (83) 



<£ kl 7t k2 £ k3 > = - (2vr) 4 (2 7T)- 3 / 2 ^^S^ + k 2 + k 3 ) 
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In this paper, we shall restrict our attention to evaluating the bi-spectrum in the 
Starobinsky model in the equilateral limit (i.e. when kx = k 2 = k 3 ). In such a case, the 
above expression for / NL simplifies to 

r = -^(2 7r)- 4 ^^ (84) 
with G(k) given by [cf. Eq (72)] 

7 

c=i 

6 

= K [fk foe) Gcik) + f* k 3 (ifo) ^(fc)] + G 7 (*). (85) 
c=i 

Using this expression and the dimensions of the quantities involved, which we had 
pointed out before, it is straightforward to check that / NL is indeed a dimensionless 
quantity. Further, in the Starobinsky model, due to the difference in the dynamics before 
and after the transition, the integrals in Eqs. (73)-(78) need to carried out separately 
on either side of the transition. Therefore, we can write 

G(k) = Ml J2{f! fa) [G + C {k] + g-{k)] + if (ifc) [g+*(k) + G-*{k)] } 

c=i 
+ G 7 (k) 

= Ml J2 [fl S+(k) + ff (ife) Q?{k)] 
c=i 

6 

+ mi [ft fa) s-(k) + /f (ifc) g-*(k)] + G 7 (k) 

C=l 

6 

= 5][G+(A;) + G C (A;)]+G 7 (A;), (86) 
c=i 

where, as we had mentioned earlier, the plus and the minus signs in the super-scripts 
refer to the quantities before and after the transition, respectively. In the following 
section, we shall first evaluate the contribution to the bi-spectrum due to the supposedly 
dominant term G4 in the Starobinsky model. And, in the subsequent section, we 
shall evaluate all the remaining contributions as well. As we have pointed out before, 
remarkably, we are able to evaluate the bi-spectrum in the equilateral limit completely 
analytically without any further assumptions or approximations. We shall then discuss 
the range of values for the parameters of the Starobinsky model for which the non- 
Gaussianity parameter / NL can be large. 

4. The dominant contribution to the bi-spectrum in the Starobinsky model 

In this section, we shall evaluate the contribution due to the term G<± to the bi-spectrum 
in the Starobinsky model in the equilateral limit. 
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To begin with, recall that, in the Starobinsky model, the second roll parameter e 2 
is constant before the field reaches the discontinuity at 0o when t] = i] = — (ao Hq)^ 1 = 
k^ 1 . As a result, the integral Q<± which involves e' 2 [cf. Eq. (76)] vanishes before r/o and, 
hence, we are left with only the following contribution after the transition: 

&(ki, k 2 , k 3 ) = i / dr a 2 ei „e 2 _ /£ + two permutations) . (87) 

J-k- 1 

Also, post-transition, the mode v k and the quantity z are given by Eqs. (29) and (30), 
respectively. Since f k = Vk/z, one obtains that 



fk(v) 



iH a k 



2M PlV /F^ 



(l+ik V )e- ikv 



iH /3 k 



2M PlV /F 



[l-ikr))e ikr >, (88) 



while the corresponding derivative with respect to the conformal time coordinate is 
given by 



iH a k 



2M PlV /P^ 
iH /3 k 

2M PlV /A?V 



-He!- (1 + ikrj) 



— Hei- (1 — ik-q) — 



ne 2 . 



[1 + ikrf) + k 2 r] 



[1 — ikrj) + k 2 T] 



-ikrf 



ikrf 



(89) 



In this expression, one can ignore the first term in the square bracket since it is 
proportional to the first slow roll parameter which, as demonstrated before, always 
remains small even when the field is crossing the discontinuity in the slope of the 
potential. However, we should stress here that, in doing so, we are in fact ignoring 
contributions to G4 that are possibly of the same order as the sub-dominant terms, such 
as Gi, G2, G3 and G-j. We shall comment more on this point in the next section. 

We find that, upon using the expressions (14) and (19) for E\ and e 2 after the 
transition, in the equilateral limit, we can write Q4 as follows: 



81 AA kl Hi 



a* k 6 h(k) - P k 12(h) - a* k 2 (3* k J 4 (k) 



2 V2¥ M 2 A 2 _ 

+ a* k pfj;(k) 

The quantities I4 and J4 in the above expression are described by the integrals 



(90) 



h(k) 



r° drr 
(1 -p 3 r 



l -ikrf k 2 + 3p 3 (1 -Oct) t 



-k 2 p 3 r 3 




c 

dr t 



3ikr 



'l-ikr) 3k 2 + 9p :i T + ik 3 T + Qp 3 k 2 T 3 



(1-P 3 r 3 ^ 



-tk 3 p 3 r 4 



Sk'i 



where the quantity p is defined by the following expression 
o 3 = -—k 3 



(91) 

3 

(92) 
(93) 
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Figure 4. The three poles corresponding to the factor (l — p 3 t 3 ) that appears in 
the integrals I4 and J4 [cf. Eqs. (91) and (92)] are marked with black asterisks in the 
complex r-plane for the case wherein AA > 0. Evidently, they will all fall on a circle 
of radius = (A-/\ AA\) 1 which is centered at the origin [cf. Eq. (93)]. Note 

that, in the integrals I4 and J4, the path of integration runs from — k^ 1 to zero along 
the negative r-axis. Amongst the three poles, viz. n, T2 and T3 [in this context, see 
Eq. (143) below], only T2 lies on the axis. The dashed blue curve represents a circle 
of radius k^ 1 about the origin. However, since |AA|/A_ < 1 when AA > 0, all the 
poles, including T2, lie outside the blue curve. Therefore, no pole falls on the path of 
integration. 



A concern could arise that, in the above integrals, one may encounter poles 
corresponding to the factor (1 — p 3 r 3 ) 1 along the path of integration, viz. from — k^ 1 
to zero on the negative r-axis. It is easy to establish that this does not occur. Since 
r is negative over the region of interest, when AA < 0, p 3 > 0, so that p 3 r 3 < 0. So, 
(1 — p 3 r 3 ) > for all r < and, hence, no pole occurs on the negative r axis in such a 
case. [The poles when AA < are given by Eq. (144), and their actual locations in the 
complex r-plane can be found represented in Fig. 7.] Whereas, when AA > 0, p 3 < 0, 
and p 3 r 3 > for r < 0. It is then clear that, under such conditions, a pole will indeed 
arise on the negative r-axis. However, as we have illustrated in Fig. 4, the pole on the 
axis of integration only occurs at r < — k^ 1 , which lies beyond the lower limit of the 
integrals. Therefore, we encounter no poles on the path of integration. 

As no poles arise, we find that the integrals I4 and J 4 can be evaluated easily, and 
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h(k) 



1 1 

3 ~ 3 



1 + i k/ko 
1 + {Plhf 



3 — 3 ik/kg 



= 1 _ (1-ik/kp) {l+ik/k ) 2 e _. fc/fco 



27 

(94) 
(95) 



[1 + {p/hff 

Upon substituting the above results for the integrals J 4 and J 4 in Eq. (90), we obtain 



that 



81 AA /cq Hi I a n 



* 3 



af f3* k e~ 2ik/ko 



3 [1 + (p/k y 



5*3 



3 



(96) 



where, for convenience, we have introduced the quantity 

P k = f3 k e- 2lk ' h \ (97) 

As we shall see, working in terms of (3k allows us to isolate exponential factors easily. 

We next need to evaluate the quantity G±{k). According to our prior discussions 
[see in particular Eqs. (86)], this quantity is given by 

G 4 (k) = Ml [f k ( Vc ) G 4 (k) + /* 3 fa) Ql{k)] . (98) 

Towards the end of inflation, i.e. as r\ — » 0, the mode reduces to 

iH 



fk(Vc) = 



(99) 



where ei_(?7e) denotes the value of the first slow roll parameter at late times. On using 
this expression and the relation (96) for C/ 4 , we find that we can write G 4 as 



k 6 G 4 (k) = 



81 



A; \ 3 AAHp 



8^2eJJri\kJ AlM^l 
+ A2 (k) cos 



Ai(k) sin 



f3k 

As(k) sin I — ) + Ai{k) cos 



V ^0 



3A- 



(100) 



where the four coefficients Ai, A2, A3 and A4 are found to be 



Ai(k) 



3A 3 _ AA 
AA% 



A + [ -9 - 



A 2 (k) = - 



3A 3 _ AA 
AA\ 



k 

9_P 
/i 



1 + 



2\ 2 



9A 
2 A; 4 

K 



1 + 
k_ 

k 



(101) 



2\ 2 



7,2 

A 



9 A_ 1 + 



/V 2 



9 + 



/t 



-43(fc) = - 



A 3 
12 



/c 2 



(102) 



/V 2 



27 AA 1 - 72 - 27AA (5A_ - 7 A + ) 



kp 
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- (9A_ - 1L4 + ) 2 ^ + QA + (-3A_ + hA + ) ^ 



(103) 



A A (k) 



A 3 _ 
\2A\ 



- 27 - 4 -(- 3+ S) ( i+ l) 2+i8 - 4 - 4+ ( i+ S 

7 k 2 6k 4 \ / 153 A; 2 9 A; 4 93 A; 6 



X r 9 "^ + ^ 4 "j + A+ 81+ A; 2 A; 4 A;^ 
4k 8 

+ 



US 



(104) 



k 

k , 

Let us now understand the behavior of at large and small scales. Experience with 
the slow roll results suggest that, far from k , on either side, one can expect A; 6 G 4 to 
turn scale invariant. As k/k — > 0, we find that 

l im k ^ G4{k) = __ J! A jf- H ° (105) 

and, in the limit A;/A; — > oo, one arrives at the result 

lim A; 6 G 4 AO = ; ° — cos — J . 106 

In Fig 5, we have plotted the absolute value of the quantity k 6 G 4 given by the 
expression (100). Clearly, while the quantity is strictly scale invariant for k <^ k , 
the scale invariant amplitude is modulated by superimposed oscillations for k 3> A; . 

It is interesting to note that, while the power spectrum had depended on 
trigonometric functions involving 2k/ko [cf. Eq. (39)], the bi-spectrum depends on 
trigonometric functions involving 3k/ko. Definitely, these behavior can be attributed 
to the fact that the power spectrum and the bi-spectrum depend on the second and the 
third powers of the curvature perturbation, respectively. In fact, that such a behavior 
can broadly be expected to occur has been pointed out previously in the literature while 
investigating a model that, in some aspects, is similar to the Starobinsky model [28]. 
The model considered earlier contains a discontinuity in the potential itself rather than 
in the derivative of the potential. Despite this difference, the slow roll parameters seem 
to behave just as in the Starobinsky model, with the first slow roll parameter remaining 
small through the transition as the field crosses the discontinuity, whereas and e'2 
grow large for a short period around the time of the transition. Since the background 
behavior is rather similar, it is possible for the two cases to be meaningfully compared. 
We find that, in order to evaluate the / NL in the model, the earlier work [28] simply 
represents e' 2 as a Dirac delta function at the transition which, then, permits an easy 
integration of the dominant contribution. However, such an approach cannot seem to 
reproduce the detailed pattern that we have obtained in our calculation here, and one 
arrives at only a rough oscillatory behavior on the small scales (in this context, see 
Eq. (6.32) of Ref. [28]). Moreover, the method does not seem to allow the calculation 
of the other, sub-dominant, contributions, as we are able to carry out in the case of the 
Starobinsky model (see our discussion in the following section). With regards to the 




Figure 5. The absolute value of the quantity k e G4, as given by Eq. (100), has been 
plotted as a function of k/ko (the blue curve). We have worked with the same of 
values of A+, A- and Hq as in the earlier figures. We should stress the fact that, as 
in the case of the power spectrum, the quantity k 6 G4 depends on the wavenumber 
only through the ratio k/ko. The green and the red curves in the inset represent the 
asymptotic behavior for k -C kg and k fco, respectively. 



overall amplitude of the effect, it is not straightforward to compare the results, since 
the models are not exactly similar. However, the order of magnitude of / NL is found to 
be in agreement with what we obtain (see Sec. 6). This leads us to conclude that our 
results are broadly consistent with those of Ref. [28]. 

5. The sub-dominant contributions to the bi-spectrum 

In this section, we shall arrive at analytic expressions for the other, sub-dominant 
contributions to the bi-spectrum in the equilateral limit. Unlike the dominant term G4 
wherein no contribution arose before the transition (due to the vanishing 62), the other 
five terms involving integrals — viz. G\, G%, G3, G5 and Gq — contribute before as well as 
after the transition. As one would expect, the calculation of the contribution until the 
transition to fast roll largely follows the computations in the slow roll case. In contrast, 
the computation of the contribution post-transition proves to be more involved, but, 
as we shall see, tractable. The last term G 7 is, obviously, straightforward to evaluate 
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since it does not involve any integrals, and only requires the amplitude of the curvature 
perturbation at late times. 



5.1. The contribution due to the second term 

We shall first evaluate the contribution due to the second term in the interaction 
Hamiltonian (67), viz. the quantity G 2 arising due to the integral Q 2 [cf. Eq. (74)], 
as it involves simpler integrals. 



x 



5.1.1. Before the transition Since e\ remains a constant before the transition, the 
integral (74) can be written as 

£?2"(ki, 1*2, k 3 ) = — 2i (ki • k 2 + two permutations) e\ + 

[ k ° dr a 2 f* kl fi 2 f* k3 , (107) 

J — (1— i 7) 00 

where the 7 which appears in the lower limit is an infinitesimal positive quantity that 
has been introduced by hand, as is done in the standard slow roll case [24, 27]. The 
procedure can also be viewed as though the integration contour has been rotated by a 
small angle 7 in the complex r plane. The inclusion of 7 ensures the selection of the 
correct choice for the perturbative vacuum. Also, algebraically, its introduction allows 
us to avoid the increasingly rapid oscillations at very early times (i.e. as rj — > —00), as 
it acts as an exponential cut-off that leads to the rapid convergence of the integral. 

Note that, before the transition, the mode v\. and the quantity z are given by 
Eqs. (27) and (30), so that we have 

f+(ri) = lH ° (l + ikri) e~ ik \ (108) 

Upon substituting this mode in the above expression for Q2, one finds that the resulting 
integral can the easily evaluated even in the most generic case of k x 7^ k 2 7^ k 3 . We 
obtain the contribution until the transition to be 

S 2 + (ki, k 2 , k 3 ) = ^J /u Tu ,s ( k i • k 2 + k, ■ k 3 + k 2 ■ k 3 ) e-^o 



4M p 3y (h k 2 k 3) 

x 

fc T Kq k T k T 

(109) 



' fcifc 2 fc 3 .ihhh.i,,, . , , 1 7 7 \ 
_ I I H V2 h T~ 2 + fc 2 % + h k 3 ) 



where k T = A7 + k 2 + k 3 . It can be easily checked that this expression reduces to the 
standard slow roll result as ko — > 00 [24, 27], which is basically equivalent to assuming 
that the transition to fast roll does not take place at all. On restricting to the equilateral 
case, the above expression simplifies to 

3H k 2 k fe^ ( 10 i k k 2 
AMI V"^ V + 9k o ~3fco, 



S}(k) - ^5 sl^ ( 1 + ^ - 4u ) e-^° (HO) 
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With the help of the above Q% , we are now in a position to calculate G%, which is 
given by 

Gt(k) = Ml [f k (rk) Qt{k) + ff (rk) g; + {k)] 

= — ~'f°3 7 s [( a k-Pk) 3 GZ(k) - K-(3l) 3 Gt(k)} ■ (Hi) 

Upon using the expressions (37) and (38) for a k and straightforward manipulations 
allow us to write G\ as follows: 



k 6 Gt(k) 



1 



ko A+H* 



16 v^U%) fc M p. 
x |3K 
10 i k 



+ 



9k 
- 33 
10 i k 



„2 o _ * 5*2 



sin 



+ 



9k 



(a 2 J k + a*Jt 2 ) (l-^|) 

„2 o _ * 5*2 



- 

+ 3 



(« 3 fc + ft*) (l - 



3 &o 



+ 



10 i k 
~9~k~o~ 



^ 3 ) 



3k 2 . 



+ 



lOiife 
9^T 



(«2 - ^ 3 
*2 - ^ 3 ) 



sin 



cos 



3A; 

k 
3 k 
k 



(112) 



where 3ft and 3 denote the real and the imaginary parts of the arguments. Notice that 
the above expression contains trigonometric functions of arguments k/k and 3 k/ko, as 
the contribution due to the fourth term had. Therefore, this expression can be expected 
to lead to a similar (but, in some ways, different!) oscillatory pattern as the dominant 
contribution. 

Let us now understand the asymptotic forms of k^G^- As k/ko — > 0, we find that 
it goes to a constant given by 

17 



lim k Gt(k) 

k/k ->0 



A 3 _ Hi 



144 y/2 4_(ri c ) M^A\- 
While, in the limit k/ko ~^ °°> one obtains the following behavior 



(H3) 



lim k 6 G+(k) = 1 

fc/fco->oo 16 y 2 ef_ (r) e 



A + H% 



1 



A. 



cos 



1 k 

+ - — sm 
3 k 



3k 
k 



47 3 A 
18 



cos 



3k 
k 



(114) 



It should be pointed out that the second term grows with k as (k/ko) sin (3 k/ko)- 111 
other words, k 6 G\ diverges linearly at large k. Interestingly, as we shall illustrate, this 
diverging term will be canceled by a similar term but with the opposite sign that arises 
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due to the contribution post-transition, ensuring that the complete k & G2 remains finite 
at large k. 

5.1.2. After the transition Let us now turn to the evaluation of the contribution 
Q2 after the transition. The calculation proceeds in a manner very similar to the 
computation of Q± . But, in contrast to Q±, for evaluating Q 2 , we only require the 
behavior of the first slow roll parameter after the transition, which is given by Eq. (14). 
Upon substituting the modes (88) after the transition in Eq. (74) and, after some 
manipulations, one obtains that 

+ alPfr 2 (k)], (115) 
where I2 and J2 are described by the integrals 

«^viife/;,S (1 - A3)(1 - i ^ )V "' (116) 
Mk) = vm^, £, % (1 - " 3t3) (1 - lkT?{1 +lkT) ^ (117) 

with p 3 given by Eq. (93). It is clear that the structure of these integrals is very similar 
to that of I4 and J4 [cf. Eqs. (91) and (92)]. However, in contrast, the poles in I2 and 
J2 are located at the origin and, hence, they do not fall on the integration path. Due 
to this reason, the two integrals I2 and J2 can be performed easily, and the results can 
be written in the following form: 

h(k) = J^-*2 Mvx [l a 2(k)+I b 2 (k)e- 3lk/k0 ] , (US) 

Mk) = * A m k l [J2(k) + J 2 b (k) e~ ik ^} , (119) 

where the scale dependent quantities I 2 , I 2 , J 2 and J 2 are given by 

^—^ski 90 ^ 53 ^ 1 )- (120) 

m = _ 1 _J_( 90i+ 53A ' 4 ^ 53iA - 4fc » 



81k \ A_k 3 J 27 A_ k 

k 2 ( AA\ 13iAAk 22AA . . 

+ 3*i( 1 -ir) + ^o^ + Tir- (121) 

Ji(k)=- 1 + ^(2, + 27AAk3A ^ 27 < AAk ° 



/cq V yi_ k^ J /uq y4._ k 

14 AA 5iAAk AAk 2 



A A. ho ' A /,vT (123) 
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Note that the coefficients T 2 and J 2 actually diverge in the limit r) e — » 0. 
This implies that I 2 and J 2 will diverge too, but, as we shall soon demonstrate, the 
corresponding contribution to Q 2 remains perfectly finite. It is also worth remarking 
that the coefficients I 2 , I 2 , J 2 and J 2 depend on the wavenumber only through the 
ratio k/ko, as expected. Using the above expressions, one obtains that 



g 2 -(k) 



A_ k 2 kn 



lZ(k)a* k 3 -3J 2 a (k)»l 2 P* k e- 2ik / k ° 

l^*{k)(3* k 3 e- 6ik/ko 



+ JC 2 {k)e- 3ik/ko + 3J 2 a *{k)a* k (3* k 2 e~ 4lk/k " 



(124) 



where the new coefficient K, 2 is given by 

}C 2 (k)=l b 2 (k)al 3 -l b 2 *(k)P* k 3 + 3a*Jl [j 2 b *(k)P* k -J 2 b (k)al\. (125) 

We can now evaluate the resulting G 2 ~ . The remarkable fact is that, all the terms 
proportional to l/i] c cancel out and, therefore, the final result for G 2 ~ turns out to be 
finite in the limit r] e — > 0. We find that G 2 can be expressed as 



k 6 G 2 {k) 
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(126) 

Let us now study the asymptotic behavior of G 2 . In the limit k/ko ~^ 0, we find 
that k 6 G 2 goes to a constant given by 

3 



A 5 FT 2 

lim k°G 2 (k) = 



(127) 
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Whereas, as k/k — > oo, it has the following behavior 
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As we had mentioned earlier, k 6 G 2 too diverges at large k. Upon comparing with 
Eq. (114), it is clear that the divergence involving the sin(3/c//co) term in the above 
expression is of the same form, but with an opposite sign. Therefore, the total 
contribution due to the second term, viz. G 2 (k) = G 2 (k) + G 2 (k) remains finite even 
at large k. We find that, as k/ko — > 0, k & G 2 turns strictly scale invariant with the 
amplitude 

lim k 6 G 2 (k) ' 



A 3 _ Hi 



17 3A 2 _ 



(129) 
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In the limit of large k/ko, the diverging terms cancel exactly and we are left with the 
following finite expression: 
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(130) 



which is scale invariant with superimposed oscillations. Notice that the superimposed 
oscillations again involve a cosine function with the argument 3 k/k . These asymptotic 
behavior are also evident from Fig. 6, wherein we have plotted the absolute value of the 
quantity k & G 2 - 



5.2. The contributions due to the first and the third terms 

Barring the k dependent factors, the contributions Q\ and due to the first and 
the third terms in the interaction Hamiltonian (67) involve similar integrals [cf. (73) 
and (75)]. One finds that, their contributions can be combined and evaluated together 
in the equilateral limit. As in the case of the second term, we shall first evaluate the 
contribution before the transition and then turn our attention to the calculation of the 
contribution after the transition. 



5.2.1. Before the transition To begin with, note that, before the transition, the mode 
fk is given by Eq. (108), while its derivative can be computed to be 

As is usually done while working in the slow roll approximation, we shall ignore the 
first term within the square brackets involving t\ + and e 2 +. Upon using the mode (108) 
and its derivative above in the integrals (73) and (75), and upon modifying the lower 




Figure 6. The absolute value of the quantity k e G2 has been plotted as a function of 
k/ko (the blue curve). In plotting this figure, we have worked with the same values 
of the parameters that we had considered in the previous figures. As in the last plot, 
the green and the red curves in the inset represent the limiting forms for k <C fco and 
k 3> ko, respectively [see Eqs. (129) and (130)]. 



limit of the integrals as in Eq. (107), we find that the integrals can be evaluated easily 
without actually having to take the equilateral limit. We obtain that 
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which reduce to the standard slow roll results as k Q — > 00 [24, 27]. In the equilateral 
limit, the above expressions simplify to 
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so that we have 



2iH k 3 
4 Ml 



k 9 



Qt{k)+Qt{k) 



iH Q k 3 
' 4 Ml 



4 

3 k 



k 9 



-3ik/ko 



4 ifc 
3 k 



-3 ik/ko 



36 



(135) 



(136) 



The corresponding + G% is given by Eq. (86) which involves the mode function 
at late times [cf. Eq. (99)], and thereby the Bogoliubov coefficients. Upon using the 
expressions (37) and (38) for a fc and (3k, we find that we can write 
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A + H% 



x 



48v^t(^ 
33 

4 
3 



Ml 



33? 









+ 






% k 


+ 




k 



+ 3 



+ & 



4 i fc 
3 

4 ifc 

3 A; 



0% Pk + «fc 

4 3 ) 



sin 



cos 



sin 



(«2 - 4 3 ) 



cos 



3 fc 
A; 
3& 
A; 



(137) 



We again encounter the by-now usual structure containing sine and cosine functions 
with arguments k/ko and 3 k/ko- As we had carried out earlier, it is interesting to study 
the asymptotic behavior of the above expression at small and large scales. In the limit 
of small k/ko, one has 



lim k e [Gt(k) + Gt(k)]=- 

k/ko—tO 

while in the limit of large k/ko, one obtains 
lim k 6 [Gt(k) + G+(k)] = 

k/ko— 
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It should be mentioned here that, as in the case of k 6 G\, the above asymptotic form too 
diverges linearly at large k. Again, as in the earlier case, the contribution post-transition 
will cancel this divergent term leading to an overall finite k 6 (G 1 + G 3 ). 



5.2.2. After the transition We shall now evaluate the contribution due to first and 
the third terms arising from the post-transition phase. As the calculation proves to be 
somewhat heavy and lengthy (but manageable), we shall relegate some of the details of 
the calculation to an appendix. The starting point is similar to what we have already 
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encountered in the case of the first sub-dominant term: upon substituting the mode (88) 
and its derivative (89) in the integrals (73) and (74), we find that, we can write 

g-(k) + g^(k) = 4 ^° 9/2 [a? i 13 (k) - fif r 13 (k) - <*? fit Mk) 

+ alfil 2 J* 13 (k)], (140) 
where I± 3 and J 13 are described by the integrals 

I 13 (k) = -==^= r < K~^ (1 - ikr) \3 p 3 r (1-ik r) 



iW Pl J -K 

+ k 2 (1 -p 3 r 3 )j 2 e 3ifcT , (141) 

[3A(l-**r) + *> (l-,^)] 

x [9 (1 -iifer) (l + iibr) p 3 r + /c 2 (l - p 3 r 3 ) (3 - i fc r)] e* kr . 

(142) 

In arriving at the above expressions, we have again ignored the term involving e±- within 
the square brackets in Eq. (89), just as we had done while evaluating the dominant 
contribution G4. We had earlier mentioned that, in the case of G4, the neglected terms 
can be expected to be of the same order as the contributions due to Gi, G2, G3 and G7. 
The terms that we shall ignore here are possibly of the same order as the contributions 
due to G5 and Gq. 

Let us now analyze these integrals in some detail. Firstly, a concern could arise that 
these integrals may contain a pole — due to the (1 — p 3 r 3 ) - term — along the path of 
integration, which runs along the real axis from — k^ 1 to r/ c , with the latter approaching 
zero from the negative direction. Actually, the integrals I4 and J 4 [cf. Eqs. (91) and (92)] 
that we had to carry out earlier in the calculation of G4 had also contained the same 
term, and we had discussed the issue of the poles in some detail then. In the case 
AA > 0, the three poles corresponding to (1 — p 3 r 3 ) 1 are actually located at 

qItt/3 gi7r g5i7r/3 

r i = -rr> r 2 — 77' Ts = ~rr> ( 143 ) 
\p\ \p\ \p\ 

in the complex r plane. In fact, we had illustrated the positions of these poles for a 
situation wherein AA > in Fig. 4. Whereas, when AA < 0, the poles are found to be 

at 

Y e 2i7r/3 e 4iw/3 

n = r-f, r 2 = ——, r 3 = (144) 
\P\ \P\ \P\ 

and their location is illustrated in Fig. 7. It should be clear from these two figures that 

no poles arise on the path of the integration. 

Secondly, in order to compute the integrals, it is convenient to write the term 

containing the poles as follows: 

1111 1 



1 — p 3 t 3 p 3 r 3 — 1/p 3 p 3 (r — T±) (7 — 7 2 ) (7 — 73) 
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Figure 7. The poles corresponding to the (l — p 3 t 3 ) term that appears in the 
integrals /13 and J13 [cf. Eqs. (141) and (142)]. We should mention here that, we had 
encountered the same term in the integrals I4 and J4 as well and, in Fig. 4, we had 
illustrated the positions of the poles when AA > 0. The location of the three poles 
Tx, T2 and r 3 in the complex r-plane, for a case wherein A A < 0, have been marked 
in the above figure with black asterisks [in this context, see Eq. (144)]. The solid blue 
line segments denote the three new integration paths corresponding to T\ , T2 and t 3 in 
the complex t n plane, with t n being related to r„ through Eq. (151). The solid green 
line running along the negative real axis represents the branch cut associated with the 
function E\(z). 



W ( 61 +_*!_ + _!»_), (145) 



p 3 \r - Ti t - t 2 t - r 3 



where 61, b 2 and 63 are given by 

p?(Tl - TMT — r -\^-^)> w 

AHr^-rMr^r-l^)^ <148) 
These expressions are valid when AA < 0. If AA > 0, then one has bi = 
— [(l + i\/3) /6], 62 = (1/3) and 63 = — [(l — i v3j /6], i.e. simply a permutation of 
the previous values. 

Lastly, the numerators of the integrals /13 and J13 are polynomials in r of order 



* 2 = 



> 3 
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seven. As a consequence, one can write 

= - j-j/i 3 y ^ r dr — ^ ( w9 ) 

a/18 # 2 M P1 p 3 ^ y_ fe -i r-r n 
where -Pi(t) is the seventh order polynomial given by 

P^t) = (1-ikr) [3p 3 r (1 -ikr) + k 2 (l -p 3 r 3 )] 2 . (150) 
Let us now perform the following change of variable from r to t n : 

t n = -3ik (r - T n ) . (151) 

It should stressed here that, because the change of the variable depends on r n , it is 
different for the three integrals in the above expression for J 13 . The change of the 
variable leads to 

/isW = - J^'l^ y b n e 3lkT - 



18^M PlP 3 r : 

x / , T e "" Pi Uri +T " ■ (152) 

J3ik/k +3ikr n l n V° 1 K J 

We find that this integral can be evaluated by initially expressing P 1 as a polynomial 
in t n , and then carrying out the integral term by term. Since the computation proves 
to be somewhat longwinded, as we had mentioned before, we shall relegate the details 
of the calculation to an appendix (see Appendix A.l), and simply quote the final result 
here. We obtain that 

A- \p\ 2 k 4 
^H 2 M Pl p 3 

where the scale dependent quantities Z° 3 and T\ z can be expressed as 

3 7 

Zi3( k ) = - y, bn [ Cn0 ^ ( 3 ^ r «) + y c «™ e m-i (3ikr n ) , (154) 



Uk) = - £M ; 3 [lUk) +4We- 3a/ 1 , (153) 
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v-^ r f3ik \ 
x i3{k) = y b n |c n0 Si + 3ikT n j 



n=l 

7 



+ y c nm (m - 1)! e m _! ( ^ + 3i A;r n ) , (155) 



m=l 



with the coefficients c nm being given by Eqs. (A.2)-(A.9). In these expressions, e n (z) 
denotes the exponential sum function (A. 14), while S±(z) = e z Ei(z), with £1(2) being 
the exponential integral function (A. 15). Despite the fact that the calculation is more 
involved, it is clear that the structure of I 13 resembles that of I 2 - 

However, there is one particular point that we have overlooked in the above 
considerations. When one performs the change of variable from r to t n — —3 ik (r — r n ), 
one basically modifies the contour along which the integrations are performed. Since 
the change in the variable from r to t n involves r n , the modified path evidently depends 
on r n . The effect of multiplying by 3ik is to rotate the path by an angle of 7r/2 in 
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the anti-clockwise direction. Therefore, the contours start from 3ikr n (corresponding 
to r = 0) and run vertically until 3ik/k has been added to 3ikr n . The three new 
integration paths corresponding to Ti, t 2 and 73 are displayed (as the solid blue line 
segments) in the complex t n plane in Fig. 7. Note that the function E\ (z) has a branch 
cut running from —00 to (represented as the solid green line in the figure). While the 
paths associated with t\ and t 3 are further away, there seems to be the possibility that 
the path corresponding to t 2 may cross the branch cut. It is clear from the figure that 
this will occur if the imaginary part of t 2 corresponding to r = — k$ 1 is positive, i.e. 
when 



9 [h (r 



3k 
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1 ( A. 
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1/3' 



> 0. 



(156) 



For AA < 0, this condition implies that A + /A_ > 9/8 = 1.125. If this condition is 
satisfied, then one crosses the branch cut of the function Ei(z) and, in such a case, 
one must add a suitable contribution to the integral ii 3 . We find that the additional 
contribution amounts to redefining the above X" 3 by 

3 r 7 1 

(3ikr n ) 



Z?s(*) 



n=l 



m=l 



+ {2tn)b 2 c 20 e 3ikT2 . (157) 

It is important to point out here that whether the additional term arises or not depends 
only on the ratio A + /A_ and not on the wavenumber k. 

We now need to calculate the integral J13, which, as in the case of J13, can be 
expressed as 



Jis(k) = 



A-\ P \ 



18 H 2 M P1 p 3 



n=l 



k a 



(158) 



where P2(j) is a seventh order polynomial given by 

P 2 (t) = [3p*T(l-ikT) + k 2 (l-p 3 r 3 )] 

x [9 (1-ikr) (1 + ikr) p 3 r + k 2 (l - p 3 r 3 ) (3-ikr)] . (159) 

To evaluate the above integral, we can now apply the same strategy that we have followed 
above for J 13 (for details, see Appendix A.l). We obtain that 
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while the coefficients d nm are given by Eqs. (A.20)-(A.27). Just as the structure of ii 3 
had resembled I2, evidently, the form of J 13 resembles that of J 2 . We again encounter 
here the issue regarding the path associated with t 2 crossing the branch cut of the 
exponential integral function E\(z). It is straightforward to show that the previous 
arguments apply for the present case as well. As a result, when A + /A_ > 9/8, instead 
of i7i3 above, one has 

3 r 7 1 

Jn( k ) = ~^2 b n d no £1 (ikr n ) + d nm (m - 1)! e m _i (ikr n ) 

n=l L m=l J 

+ (2i7i)b 2 d 20 e ikT2 . (163) 

The above expressions for J 13 and J 13 then determine C/f + Q 3 [cf. Eq. (140)], 
which in turn can be used to arrive at the corresponding Gj~ + G 3 . We find that the 
calculation proceeds in exactly the same fashion as in the evaluation G 2 , because, as 
we have already pointed out, the final expressions for the integrals I 13 and J 13 are very 
similar to that of I 2 and J 2- As a consequence, one can write 
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x (a k al-~h~^i)} cos (^))> ( 164 ) 



where the coefficient /C13 has been defined to be 



/Cis = l b 13 (k) a* k A - l b 13 *(k) ft? + a* k PI [J^ik) PI - JUk) a* k \ . (165) 

Obviously, in the calculation of the quantity G~[ + G 3 , the coefficient /C13 plays the 
same role that JC 2 had played in the evaluation of G^. Therefore, it does not come as 
a surprise that the structure of /C13 is identical to that of K 2 - Moreover, the expression 
of G\ + G 3 itself bears a lot of resemblance to G2 with, in particular, the 'harmonics' 
k/ko, 2 k/ko, 3 k/ko and A k/ko, all being present. 

Let us now calculate the limiting forms of the above expression for small and 
large k/k . As k/ko — > 0, we find that G± + G 3 behaves as follows: 

I A3 tt2 / u 

hm k % \G^(k) + G 3 (k)] = ( T- 

k/ko^o 3V ;J 320 v /2e?_(7 7e ) M^A 2 + \k 

( 2A_ 3A 2 \ 

x ( 15 + — + (166) 

i.e. it goes to zero quadratically for small wavenumbers. Therefore, in this limit, it is 
the term before the transition that dominates and, hence, the complete contribution 
due to the first and the third terms reduces to 

lim k 6 [Gi(k) + G 3 (k)] = lim k 6 [G+(k) + G+(k)] 

k J ko ^0 k J /so ■ ^0 



1 A 3 _ Hi 



36 ^2e(_(Ve) A\ ' 



(167) 



In the limit k/ko 00, one finds that 

lim k 6 [G^(k) + G 3 (k)] 



A- Hi 



*/*o->°o L /J 36^2 4Jjk) 



27 'AA ( k\ 36A+A; (3k 
+ — cos — — — — — —sin — 

8A_ \k J A8A^k \k 



AA ( AA\ AA 

'50 + 27— - 8 - 27 



1 

x cos ( ^ I \. (168) 



A_ \ A+J A_ 



v ^0 

As usual, only trigonometric functions with the argument 3 k/ko remain in the final 
expression of the limit. Also, note that the coefficient of the sin (3 k/ko) term diverges 
linearly at large k, exactly as the contribution before the transition had [cf. Eq. (139)], 
albeit with an opposite sign. Hence, the complete contribution remains finite and, we 
find that, in the limit k/ko — > 00 , Gi + C 3 goes to 




Figure 8. The absolute value of the quantity k 6 (G\ + G3) has been plotted as a 
function of k/ko (the blue curve). We have again worked with the same values of the 
parameters as in the earlier figures. As in the plots before, the green and the red curves 
in the inset denote the asymptotic forms. 
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(169) 



As before, at large wavenumbers, we obtain a scale invariant amplitude which is 
modulated by super-imposed oscillations of the form cos (3 k/ko). In Fig. 8, we have 
plotted the absolute value of the quantity [k 6 (G\ + G 3 )}. The figure clearly reflects the 
limiting forms that we have obtained above. 
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5. 3. The contribution due to the fifth and the sixth terms 

Let us now turn to the evaluation of the contribution due to the remaining two terms 
in the interaction Hamiltonian (67), viz. the fifth and the sixth. As in the case of Q\ 
and C/3, we find that, the terms Q 5 and Q 6 too contain the same types of integrals. Hence, 
in the equilateral limit, they can be evaluated together. 



5.3.1. Before the transition Since e\ is a constant before the transition, the integrals 
encountered in Q 5 and Q§ prove to be exactly of the same form as in the cases of Q\ 
and Qz- Therefore, they can be evaluated in the same fashion and, in the most generic 
situation of k x ^ k 2 ^ k 3 , one obtains that 
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In the equilateral limit, we find that 
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(172) 



From the last equality, we can then immediately conclude that Gf(k) + G^(k) = 
— 3e 1+ [Gf(k) + G^(k)] 1 4. As a consequence, the asymptotic behavior of G\ + G\ 
can be arrived at from the forms of G\ + G\ we had obtained earlier. As k/k — > 0, 
one has 
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AlH 2 

48v/2et(^c) M l A \ 
1 A 3 



864 V2 e?_(%) M p 7 . H2 o 



2 ' 



(173) 



lim k 6 [Gt(k) + G+(k)] 

fc/fco— >co 



1 



A 3 



1152 y/2eUVe) M l H l 
2 

3A; 



x 



1 — cos — 

A+J \k 



k 

+ jr sin 
fc V fc 



The scalar bi-spectrum in the Starobinsky model: The equilateral case 



45 



+ 



35 
~6~ 



9 A. 
21 



cos 



3k 
ko 



(174) 



Note that, while k 6 (Gf + Gq) goes to a constant at small k/ko, it diverges linearly 
(modulated by oscillations) at large k/ko, as in the case of k 6 (Gf + G^). Also, as we 
shall illustrate, just as in the earlier cases, the diverging term will be exactly canceled 
by a corresponding term that arises post-transition. 



5.3.2. After the transition The next step is to perform the calculation after the 
transition. The calculation proceeds exactly as in the case of + Q$ . We find that, 
we can write 

+ al/3l 2 J* m (k)_ . 

where I 56 and J 56 are described by the integrals 
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These integrals are very similar to the ones considered before but, crucially, they do 
not contain any poles. This makes their explicit calculation considerably easier. Upon 
following the same strategy that we had adopted in the cases of /13 and J13, one obtains 
that 
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where the coefficients f n are given by Eqs. (A.30)-(A.40), while e n (z) is the exponential 
sum function (A. 14). Similarly, we find that can be obtained to be 

j - (k) -^^y- {k)+jk{k)e '' klh ' ] (i8i) 
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(182) 
(183) 



and the coefficients g n being given by Eqs. (A.41)-(A.51). 

At this stage, the calculation of G^+Gq progresses as before and the final expression 
reads 
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185) 



where the coefficient /C 56 has been defined to be 

K x (k) = lUk) af - l b 56 *(k) fa 3 + a* k fa [j 5 b 6 *(k)fa - J 5 b 6 (k)a* kl 

Needless to add, the coefficient /C56 is similar to /C 2 and /C13. Further, as far as the 
structure of G^ + G e is concerned, it again involves the various trigonometric functions 
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that we had observed in the earlier cases. 

Let us now evaluate the asymptotic forms of the above expression. As k/k — > 0, 
we find that 
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In other words, this term vanishes quadratically for small wavenumbers, just as the 
contributions due to the first and the third terms, post-transition, had. So, in this limit, 
the total contribution due to the fifth and the sixth terms reduces to the strictly scale 
invariant value arising before the transition, i.e. we have 
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k/ko^-0 



In the limit k/k — > oo, we find that 
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If we now combine the asymptotic forms before and after the transition, for large k/k , 
we find that the terms proportional to k/k cancel, and one is led to the expression 

1 A 3 



lim k 6 [G 5 (k) + G 6 (k)} 



fc/fco^oo " - ■ - " 864 y/2^Jjk) M lx H o ' 
which is exactly the same as the scale invariant quantity that we had encountered in 
the limit of small k/k . In particular, it is interesting to note that, in this case, no 
superimposed oscillations are present. These asymptotic behavior are evident in Fig. 9, 
wherein we have plotted the absolute value of k 6 (G 5 + G 6 ). 
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5.4- The contribution due to the field redefinition 

Let us now turn to the contribution Gj due to the field redefinition (79) which in the 
equilateral limit simplifies to 
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Figure 9. The absolute value of the quantity k 6 (G^+Gq) (in blue) and the asymptotic 
behavior (in green and in red in the inset) has been plotted for the same set of 
parameters as in the previous figures. The inset highlights the fact that the quantity 
goes to the same scale invariant value at small and large wavenumbers. 



where the last equality follows from the fact that, at late times, for the linear potential 
of our interest, e 2 _ = 4ei_. Upon using the expression (99) for f k at late times, we 
obtain that 
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Using the expressions (37) and (38) for a k and (3k, we find that, we can write 
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Figure 10. The absolute value of the quantity k 6 G7 (in blue) and its asymptotic 
forms (in green and in red) have been plotted for the same set of parameters as in the 
previous figures. 



Far away from the characteristic scale k , this quantity turns strictly scale invariant, 
just as the power spectrum does. Its asymptotic forms are found to be 
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In Fig. 10, we have plotted the absolute value of the above expression for k 6 G-j. The 
figure clearly illustrates that the quantity turns scale invariant asymptotically. 

This concludes our calculation of the complete bi-spectrum. With all the 
expressions at our disposal, we can now relate these results to the observable parameter 
/ NL , and also discuss the various possible conclusions and implications. 



6. / NL in the Starobinsky model 

In this section, we shall discuss two issues. We shall firstly focus on whether the 
Starobinsky model can lead to as large a value for the non-Gaussianity parameter / NL as 
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the currently quoted mean values. Recall that, as we had mentioned in the introduction, 
the recent CMB data constrains the parameter to be / NL = 32 ± 21 in the local limit. 
However, at this stage, it is important to stress that the local non-Gaussianity parameter 
is a priori different from a scale dependent / NL in the equilateral case. In fact, the 
relevance of the constraints on the local parameter to the that we obtain below is 
not entirely clear. This point needs to be borne in mind when we compare the / NL in the 
equilateral limit with the observational constraints quoted above. Then, in the second 
part of this section, we shall focus on an issue related to the hierarchy of the various 
contributions to the bi-spectrum. 

6.1. Can be large in the Starobinsky model? 

Let us now discuss as to how large can the non-Gaussianity parameter be in the 
Starobinsky model. Since it is the fourth term that seems to often provide the dominant 
contribution to the bi-spectrum (in this context, however, see, the following sub-section), 
the quantity can be evaluated based on this contribution. On substituting the 
expression (100) for k 6 G 4 in the definition (84) of / NL in the equilateral limit, we obtain 
that 



with the power spectrum V s being given by Eq. (39), whereas the coefficients A\-A 4 
are given by Eqs. (101)— (104). Let us now consider the various qualitative aspects 
of this result || . It is important to notice that, as already remarked in the section 
devoted to the calculation of G 4 , / NL is not simply given by an expression proportional 
to sin(3/c//co + <fo), y?o being a phase (in this context, see, Ref. [28]), but contains 
four different terms that oscillate with different frequencies. It is only in the limit 
k/ko — > oo that a simple behavior as the one previously mentioned is recovered. 
Moreover, interestingly, we find that the expression for f^ 4 ^ above depends on the 
slopes A + and A_ of the potential only through the ratio R = A_/A + , and it does not 
depend on the parameter Hq at all. As should be clear by now, both the power spectrum 
and the bi-spectrum [or, equivalently, G n (k)] exhibit a step in the Starobinsky model. 
Note that, in all the earlier figures, we had worked with the parameters A + and A_ 
such that R = 0.216 < I. As is obvious from the figures, the step always rises towards 
the larger wavenumbers in such a case. In contrast, one finds that the step in the power 

|| As we were completing this manuscript, a preprint appeared, which deals with a similar but different 
model, and also evaluates the non-Gaussianities analytically [36]. We should emphasize that our effort 
is more complete, as we evaluate the complete bi-spectrum rather than just focus on the dominant 
contribution. 
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Figure 11. The non-Gaussianity parameter in the equilateral limit, i.e. f***, due 
to the dominant term in the Starobinsky model. The blue curve corresponds to the 
values of the parameters that we had considered in the previous figures. The red curve 
corresponds a larger A + (we have set A + /M^ 1 — 3.35 x 10~ 13 ), but with the remaining 
parameters being the same as for the blue curve. The dashed green lines represent the 
corresponding asymptotic values for small k/ko [given by Eq. (196)]. The insets exhibit 
the oscillations about zero, with the dashed green curves highlighting the behavior at 
large k/k . 



spectrum and the bi-spectrum reverse direction when R > 1, with the lower level of 
the step being located at large wavenumbers. The height of the step depends on the 
extent of the difference in the slopes of the potential on either side of the discontinuity. 
The further is R from unity, the larger proves to be height of the step. This property 
becomes explicit in the limit k/k — > wherein the dominant contribution to / NL due 
to the term G4 has the following simple form: 

v^» / » 4, = ^-*>- (196) 

In Fig. 11, we have plotted the expression (195) for for two different values of the 
ratio of the slopes of the potential. And, Fig. 12 contains a two-dimensional contour plot 
yeq(4) j n plane of R and k/ko- It is clear from these figures that, in the Starobinsky 
model, the non-Gaussianity parameter f** can be as large as the currently indicated 
mean value. 
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Figure 12. A two-dimensional contour plot of the non-Gaussianity parameter 
due to the dominant term in the Starobinsky model, plotted in the plane of k/ko and 
R = A-/A + . The white contours indicate regions wherein can be as large as 50. 
Note that, provided R is reasonably small, can be of the order of 25 or so, as is 
indicated by the currently observed mean value. 

But, the question that immediately springs mind is whether can be large for 
values of the parameters of the Starobinsky model for which the power spectrum proves 
to be consistent with the data. In plotting all the figures, we have ensured that the 
parameters that we have been working with lead to the COBE normalization of the 
power spectrum at suitably small scales. Given a Hq, COBE normalization restricts 
the parameter A_ to a fixed value [cf. Eq. (41)], while leaving the parameter A + or, 
equivalently, R, unconstrained. Though the actual power spectrum that arises in the 
Starobinsky model in itself has not been compared with the CMB data, comparisons of 
some variations thereof have been been carried out. These variations have essentially 
involved introducing, by hand, an overall multiplicative factor to incorporate a suitable 
spectral tilt or a sharp cut-off at the lower wavenumbers (see the fourth reference in 
Refs. [8]). While the tilt ensures a good fit to the data at the higher multipoles, a 
sharp cut-off can improve the fit to the outliers at the lower multipoles that we had 
discussed in some detail in the introductory section. One finds that, in the presence 
of a suitable spectral tilt, the best fit value of R proves to be about 0.73t° H (see the 
fourth reference in Refs. [8]; also, according to the article, the best fit value for k is 
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k ~ 3.1^2! x 10~ 4 Mpc -1 ) and, therefore, the domain R < 1 is favored by the data. 
This implies that R = 0.216 is about 3.6 ex from the mean value, whereas R = 0.00216 
is further away, being at 5.21a, from the mean value. Clearly, had these results been 
obtained for the exact case of the Starobinsky model (i.e. in the absence of the overall 
tilt that has been introduced by hand), low values of R, such as R = 0.00216, would 
be ruled out. In summary, as is evident from Fig. 11, the Starobinsky model can result 
in a power spectrum that remains reasonably consistent with the data while at the 
same time lead to an / NL that is as large as a few, at most ~ 0(10). Indeed, it would 
be interesting to strengthen such a conclusion by carrying out an explicit comparison 
of the power spectrum in the Starobinsky model with the data and also studying the 
exact implications of an oscillatory, scale dependent, / NL for the non-Gaussianities in 
the CMB. 

Though the Starobinsky model, with parameters that are consistent with the data, 
in itself, may not lead to a large / NL , we believe that the Starobinsky model with a much 
smaller R, and hence a rather large / NL can be considered to effectively capture certain 
aspects of other models that are known to lead to a better fit to the data than the 
conventional, power law, primordial spectrum. An example of such a model would be 
the punctuated inflationary scenario [10]. In punctuated inflation, a step in the power 
spectrum arises exactly as in the Starobinsky model with R < 1. However, the step 
turns out to be sharper (than in the case of the best fit value of R ~ 0.73), and the 
oscillations at the higher level of the step, before the spectrum turns scale invariant, 
are fewer. The Starobinsky model with a much smaller R can broadly be considered 
to mimic the sharper step that arises in punctuated inflation. Therefore, if one naively 
extends the results of the Starobinsky model to punctuated inflation, it suggests that 
there can exist scenarios which lead to a sharp step in the power spectrum (along with 
certain characteristic oscillations), an improved fit to the data as well as a reasonably 
large / NL . 

6.2. The hierarchy of contributions to the bi-spectrum 

As we have argued before, prior experience suggests that, when departures from the slow 
roll arise, it is the fourth term in the interaction Hamiltonian (67) that leads to the most 
significant contribution to the bi-spectrum. However, such a conclusion has largely been 
based on numerical analysis, which typically involves working with specific values (or, 
at the most, a limited range of values) for the parameters of the model concerned. Since 
we have been able to arrive at analytic expressions for all the different contributions to 
the bi-spectrum for the Starobinsky model, it is interesting to investigate whether such 
a conclusion indeed applies for a wide range of the parameters of the model or if there 
exist certain values of the parameters for which other terms can possibly contribute as 
much as or even more than the fourth term. 

Let us quickly recall some essential properties of the various contributions to the 
bi-spectrum, arising out of the interaction Hamiltonian and the field redefinition. While 
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the contributions due to the first, the second and the third terms in the interaction 
Hamiltonian involve e\ , the fourth term contains e± e' 2 , whereas the fifth and the sixth 
terms depend on e\ [cf. Eqs. (73)-(78)]. Also, note that, apart from the difference 
in the dependence on the slow roll parameter ei, the first, the third, the fifth and 
the sixth terms depend on the curvature perturbation and its derivative in the same 
fashion, viz. they depend linearly on the curvature perturbation and quadratically on 
its derivative. In contrast, the second term does not depend on the derivative of the 
curvature perturbation at all, while the fourth term depends linearly on the derivative 
and quadratically on the curvature perturbation itself. Further, the additional seventh 
term arising due to the field redefinition depends on the late time behavior of e 2 and 
the fourth power of the curvature perturbation. In a slow roll inflationary scenario, 
one can anticipate that the first, the second, the third and the seventh terms lead to 
contributions of similar magnitudes to the bi-spectrum, while the contributions due to 
the fourth, the fifth and the sixth terms can be expected to be suitably suppressed 
in amplitude due to the presence of the additional slow roll parameter. In fact, these 
expectations are broadly corroborated by the results obtained [24, 27]. 

When there arise deviations from slow roll, as we have repeatedly mentioned, it is 
the fourth term that is expected to dominate, with the first, the second, the third and 
the seventh terms proving to be of roughly similar order, but smaller in amplitude than 
the fourth term. The contributions of the fifth and the sixth terms can be expected to 
be further smaller in amplitude. These expectations are confirmed by Fig. 13, wherein 
we have plotted the absolute values of all the contributions to the bi-spectrum in the 
equilateral limit that we have focused on. It is worth pointing out here that the slightly 
higher amplitude of the second term in contrast to the contributions due to the first 
and the third terms can be attributed to the fact that the second term does not involve 
any derivative of the curvature perturbation. We shall return to this point a little later 
in our discussion. 

Since the hierarchy of the different contributions changes as one moves from the slow 
roll to a fast roll regime, it is interesting to identify the domain where the hierarchy shifts 
in the Starobinsky model. Recall that, a crucial assumption of the Starobinsky model is 
that it is the constant Vo in the potential which is dominant near the discontinuity. Once 
this condition is satisfied, the departures from slow roll can essentially be described in 
terms of the ratio R. Upon plotting the various contributions for different R, we find 
that, deviations from the slow roll hierarchy begin to occur even for an \R — 1| as small 
as 10 -5 . 

With the analytic results that we have at hand, it is also worthwhile to investigate 
whether the fourth term remains the dominant term for all values of the parameters of 
the Starobinsky model when deviations from slow roll occur or if there exist regimes 
where the hierarchy is mixed or even, possibly, absent. Motivated by such an aim, let us 
enquire if, say, the second term G 2 can probably be as large as the fourth term G 4 for 
any set or a range of parameters of the Starobinsky model. Let us consider, for instance, 
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Figure 13. The contributions due to the different terms to the bi-spectrum — viz. 
k 6 times the absolute values of G\ + G3 (in light green), G2 (in red), G4 (in blue), 
G5 + G% (in purple) and G7 (in dark green) — that we had plotted in the earlier figures 
(i.e. in Figs. 8, 6, 5, 9 and 10, respectively) for the Starobinsky model, have been 
assembled here to illustrate the hierarchy in a fast roll regime. The hierarchy of the 
various contributions to the bi-spectrum is evident from the figure. 



the ratio 



Q(k) 



G 2 (k) 



(197) 



where / eq ( 2 ) and / eq ( 4 ) are the contributions to the non-Gaussianity parameter / NL due 
to the terms G2 and G4 in the equilateral limit. As we have discussed, the first slow 
roll parameter before the transition, viz. ei+, is expected to be small. In plotting 
the earlier figures, we have worked with parameters such that e 1+ ~ 10~ 4 (cf. Fig. 1). 
But, the various assumptions and the approximations of the Starobinsky model will 
remain valid even if e i+ is larger, say, about 1CT 2 or so. If we choose ei + ~ 10~ 2 , we 
can have G2 — G4, provided Q ~ 10 2 . Interestingly, we find that the quantity Q can 
be expressed completely in terms of R and, as in the case of all the other perturbed 
quantities, its dependence on the wavenumber arises only through the ratio k/k^. In 
Fig. 14, we have plotted the absolute value of the quantity Q as a function of k/k for a 
few different values of the ratio R. It is clear from the figure that, there does exist ranges 
of parameters of the Starobinsky model for which Q ~ 10 2 and, therefore, G2 — G4, in 
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Figure 14. The absolute value of the quantity Q, which essentially reflects the 
amplitude of the contribution due to the second term Gi in contrast to the fourth 
term G4, has been plotted as a function of k/ko for a few different values of the ratio 
R = A-/A+. We should stress again that Q depends only on R and k/ko. One finds 
that, deviations from slow roll begin to occur for \R— 1| as small as 10 -5 . The upward 
spikes correspond to the wavenumbers at which G4 vanishes, so that the contribution 
due to Gi turns dominant. Even apart from these specific locations, it is evident that 
Q can be as large as 10 2 in a domain which admits fast roll. Since ei+ can be possibly 
as large as 10 -2 with the assumptions and approximations of the Starobinsky model 
continuing to remain valid, this indicates that G2 can be of the order of G4 in a fast 
roll regime. 



domains that could admit departures from slow roll (i.e. when R deviates sufficiently 
from unity). To our knowledge, this conclusion is new, since, in the literature, G4 has 
always been considered to be the dominant term when deviations from slow roll occur. 
Moreover, we believe that this is the first time that this phenomenon has been explicitly 
illustrated. This suggests that, it is plausible that the hierarchy is different in different 
regions of the parameter space of a model. However, a cautionary remark needs to be 
added to arriving at such a conclusion. Recall that, in our evaluation of £4, we had 
neglected certain terms in the quantity f' k ~ [see our remarks that immediately follow 
Eq. (89)]. As a result, a concern could arise that, for the parameters discussed above, 
the neglected terms contribute sufficiently to restore the hierarchy wherein G4 is the 
dominant term during fast roll. But, we have also analyzed the issue using a numerical 
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code [37] and our preliminary results seem to confirm the above conclusion, viz. that 
the neglected terms are not significant enough to restore the commonly expected fast 
roll hierarchy. 

In order to reinforce and also to explicitly illustrate the above conclusion, it is 
important to arrive at a specific set of parameters of the Starobinsky model which allow 
fast roll as well as lead to G 2 and G 4 that are of the same order. But, in order to 
converge on a viable set of parameters, a few points concerning certain aspects of the 
Starobinsky model and the results that we have obtained requires some emphasis. In 
particular, as k/k — > 0, we find that Q has the following simple form: 

1 n 97 R 2 

«*> = ^wk- < 198 > 

In this limit, Q ~ 17/27 ~ 0.63 for small R, while Q ~ R, when R is large. These 
behavior are clearly reflected in Fig. 14. So, if we require Q to be large, we can either 
work with a high enough R and focus on small wavenumbers, or choose an R that is 
suitably smaller than unity which permits sufficient deviation from slow roll and also 
possibly leads to G2 — G4 at large wavenumbers. However, recall that the regime R < 1 
seems to be favored by the CMB data. With these points in mind, let us now focus on 
arriving at a set of parameters in the R < 1 domain which leads to Q ~ 10 2 . Since, given 
a H , A_ is already determined by COBE normalization, choosing an ei+ fixes A + . We 
find that, if we work with the earlier values of H and A_ (listed in the caption of Fig. 1), 
while set e 1+ ~ 1/75 ~ 0.013, corresponding to A+/M^ = 3.87 x 10~ 13 and R = 0.0188, 
the hierarchy of the contributions is indeed different, with G2 being of the same order 
as G4 at large wavenumbers. This point is evident from Fig. 15 wherein we have plotted 
the various contributions G n for these values of the parameters. Of course, such a low 
value of R lies about 5 a away from the best fit value of R ~ 0.73 that we had mentioned 
earlier. Given that the Starobinsky model in itself has not been directly constrained 
(only certain variations around it have been tested precisely), this could indicate that the 
above choice of R is possibly not completely ruled out by the data. However, in the case 
of the Starobinsky model, it seems to us that, an altered hierarchy of the contributions 
to the bi-spectrum is systematically associated with regions of the parameter space that 
is not favored by the CMB data. If this is indeed true, this obviously tones down the 
importance of this modification to the hierarchy in the Starobinsky model. But, as we 
had discussed in the previous sub-section, the Starobinsky model with a smaller R may 
be considered to mimic other models that provide a good fit to the data. It is possible 
that one may encounter such altered hierarchies in these models. 

Before we conclude, let us briefly touch upon the possible reasons for the variations 
to the hierarchy. We find that the understanding gained on the hierarchy of the different 
contributions to the bi-spectrum has been largely based on the behavior of the slow roll 
parameters. Needless to mention, in addition to the slow roll parameters, the various 
contributions involve the curvature perturbation and its derivative as well. In a slow roll 
inflationary scenario, the amplitude of the curvature perturbation evolves monotonically 
once it leaves the Hubble radius. In particular, while the amplitude of the mode fk 
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Figure 15. The power spectrum (on top) and the contributions due to the different 
terms to the bi-spectrum (beiow) for a case wherein Q ~ 10 2 . The power spectrum and 
the different contributions to the bi-spectrum have been plotted in the same fashion as 
in Figs. 3 and 13. In arriving at these plots, we have worked with the same values of H a 
and A_ as in Figs. 3 and 13, but we have set A + /Afp, = 3.87 x 10~ 13 , which corresponds 
to R = 0.0188. The match between the analytical and the numerical results in the 
case of the power spectrum is rather good, which indicates that the assumptions and 
approximations of the Starobinsky model are quite valid for the parameters that we are 
working with here. Interestingly, in contrast to Fig. 13, where G4 was the dominant 
term, here we find that the contribution due to G2 is of the same order as G4 at large 
wavelengths. 
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quickly goes to a constant value, its time derivative f' k dies down as e~ 2N at super 
Hubble scales. Because of such a monotonic behavior of the curvature perturbation 
during slow roll, in such situations, the hierarchy is largely dependent on the slow 
parameters themselves. However, when departures from slow arise, it is known that the 
modes which leave the Hubble radius just before or during the periods of deviations from 
slow roll can evolve strongly around the time of Hubble exit. In fact, at super Hubble 
scales, the amplitude of these modes can be enhanced or suppressed compared to their 
amplitudes at Hubble exit (in this context, see, for instance, Refs. [30]). Therefore, 
it is seems that, when departures from slow roll occur, the hierarchy can be different 
based on a mixture of the non-trivial evolution of the slow roll parameters and the 
curvature perturbation. We believe that this is an interesting aspect that demands 
further investigation. 

7. Summary and outlook 

As we had discussed in some detail in the introductory section, models that lead to 
features in the primordial spectrum gain importance due to the fact that certain features 
can lead to a better fit to the CMB data than the conventional, nearly scale invariant 
scalar power spectrum, as is generated by slow roll inflation. The generation of features 
in the primordial spectrum requires one or more periods of of fast roll. While the scalar 
power spectrum and the bi-spectrum can be evaluated in a model independent fashion 
in the slow roll approximation, such a model independent approach seems difficult when 
departures from slow roll arise. This is essentially due to the fact that the deviations from 
slow roll can occur in a multitude of forms and it is impossible to describe all possible 
deviations in terms of a limited number of variables or parameters (for a broader effort, 
that attempts to capture a certain class of departures from slow roll, see Ref. [18]). 
In fact, in the literature, numerical computations are often resorted to in order to 
investigate scenarios involving fast roll. 

In such a situation, the Starobinsky model provides a welcome relief as it allows 
the background as well as the scalar power spectrum to be evaluated analytically 
to a very good approximation, even though it contains departures from slow roll. 
In this work, we have shown that the scalar bi-spectrum can also be computed 
analytically. Simultaneously, we have also been working towards numerically computing 
the scalar bi-spectrum and the non-Gaussianity parameter / NL in the Starobinsky model. 
Preliminary investigations indicate that the analytic expressions we have obtained match 
the numerical results quite well [37], which reflect the extent of the accuracy of the 
assumptions and approximations of the Starobinsky model. 

Interestingly, we have found that, in the Starobinsky model, for certain values of 
the parameters, the non-Gaussianity parameter / NL in the equilateral limit can be as 
large as indicated by the currently observed mean values, with the possible limitation 
that this occurs in a region of the parameter space that is not favored by the current 
CMB data. We had focused here on the equilateral limit and, clearly, it is imperative 
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that we extend the analysis to other configurations such as, say, the squeezed limit [38] . 
Further, there exist a few other models in the literature which allow the background 
and the perturbations to be evaluated analytically, and it is worthwhile to evaluate the 
bi-spectrum in these models too. We are currently investigating these issues. 
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Appendix A. Evaluation of the integrals 

This appendix contains some additional information pertaining to the evaluation of 
certain non-trivial integrals we had encountered in the text. 

Appendix A.l. Evaluation of I i3 and J 13 

In this sub-section, we shall provide a few details regarding the evaluation of the integrals 
Ji3 and J 13 . 

Let us start with the expression (141) for /13. As we had mentioned, the integral 
can be evaluated by writing P\ as a polynomial in the variable t n , and computing the 
integral term by term. If we write 

Pl(^+T n )=fc 4 5> nm C (A.l) 

we find that the coefficients c nm are given by 

c n0 = - 27 is5 3 -27e5 2 e ie " + 9 (ie5 + 5 4 ) e 2if \ (A.2) 
c nl = -42ie5 3 + (-27e5 2 + 6i5 5 ) e 10 " + (6ie 5 + 27 5 A ) e 2te " , (A.3) 
c n2 = -22ie5 3 -5 6 -9 (e5 2 -i5 5 ) e ze " + (ie5 + 22 5 4 ) e 2if \ (A.4) 

c n3 = - ^ is 5 3 - 5 6 + (-e5 2 + f l 5^ e*» + ^ 4 e 2 *, (A.5) 

c n , = - I is 5 3 - f 7 5« + 1 + § 5* e 2 *^, (A.6) 

c n5 = -^ 6 + ^ 5 e^ + ^ 4 e 2 ^, (A.7) 

c n6 = -4* 6 + ^ 5 e- (A.8) 

*r = -^, (A.9) 

where 5 = (\p\/k) = (\AA\/A_) 1/3 (k /k) and e = 1 if AA < and e = -1 if AA > 0. 
(The quantity e that we have introduced here should not be confused with the slow roll 
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parameters.) Moreover, when AA > 0, then we have Q\ = n/3, 9 2 = it, and 9 3 = 5 7r/3. 
On the other hand, if AA < 0, then we find that B x = 0, 6 2 = 2n/3 and 6 3 = 4vr/3. 
Upon using the expression (A.l), 1 13 can be written as 

12 la 3 7 



/i 3 (fc) = - 



A_ \p\ 2 k 4 



dt n C _1 e-*" 



18 H2M Pl p 3 ^ ^ n 

U "r ra =l m=0 
/>— 3ifcr7e+3ifcT, 

X / ""n f n 

•/ 3 i fc/fco+3 j fe r n 

which can be integrated to yield 

A_|p| 2 A; 4 7 



3 ?■ TV) 



(A.10) 



/is(A0 = 



18 # 2 M P1 p 3 



EE*. 



3 2 fe t"t^ 



n=l m=0 

T (to, —3ikr] e + 3ikT n ) 



v i 3 i k q ■ 7 

1 m, — h 3ikr n 

k 



where r(m, z) is the incomplete Gamma function (see, for example, Ref. [39]) 



r(m, z) 



dt e~ l t 



t j.m—1 



(A.11) 



(A.12) 



If m in a non-vanishing integer, the incomplete Gamma function reduces to an 
elementary function [39], viz. 

r(m,*) = (m-l)!e-*e m _i(;s), (A.13) 
where e n (z) is the exponential sum function defined as 



z n {z) = J2^- 



£=0 



However, if m = 0, then 

r(o, z) = e x (z) = J C 



dt 
T 



(AAA) 



(A.15) 



where Ei(z) is the exponential integral function [39]. For this reason, it is convenient 
to split the sum over to into two parts as follows: 

1 2 lA 3 



hs(k) = - 



A_ \p\ 2 k 4 



5> 



.,3 i /C Tyj, 



18ff 2 M P1 p 3 ^ 



CnO -^1 



3 i 



+ 3ikr n 



— CnO 

E l {3ikr n )+ 

"n 



171=1 



3ik 

r to, — h 3 1 kr n 

k 



T (to, 3ikr n ) 



(A.16) 



where we have taken the limit 7/ e — ?■ 0. In terms of the exponential sum function e n (z) 
and the exponential integral function Ei(z), we can rewrite the above expression for /13 

as 

12 hA r 3 



hs(k) = - 



A_ \p\ 2 k 4 



18 Hi M Pl p 3 



-5> 



c n0 e SlkT -E 1 (3ikr n ) 
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m=l 



+ e 



-3ik/k, 



n=l 



c n0 e 3ik/ko+3ikTn E 1 



3ik 



3ik 
k n 



+ 3 % k r r 



El \ I OIK 

c nm (m- l)!e m _i I + 3?fcr, 



m=l 



(A.17) 



which is the expression (153) we have quoted in the text. 

Let us now turn to the evaluation of J13 described by the integral (158). In this 
case, we can perform the change of variable t n = —i k (r — r n ), which leads to 

12 3 



Mk) = - 



A- \ P \< 



x 



18 H 2 M P1 p3 ^ 

-ikric+ikr-n ^ 

7 



n=l 

n e" 4 " P 5 



' j fc/fc() +i fc Tn 

We can write, as in the case of the function P 1 , 



i k 



^ ' m=0 



where the coefficients d nm 's are given by 

d n0 = -27^5 3 + 27e5 2 e ie " +27 (-ie5 + 5 4 ) e 2ie " 



(A.18) 



(A.19) 



(A.20) 



d n i =5Aie5 3 + (99e5 2 + 5AiS 5 ) e l6n + (-IBieS + 81 5 4 ) e 2l6n , (A.21) 

d n2 = 162ie5 3 - 275 6 + (45e5 2 + 81i5 5 ) e l6n + (9ie5 - 545 4 ) e 2ie \ (A.22) 

d n3 = 60ie 5 3 - 27 5 6 + (-27 e 5 2 - 36i5 5 ) e l6 " - 1505 4 e 2 * e ", (A.23) 

d n4 = -33ie5 3 + 9 5 6 -75 i5 5 e ie " -45 <5 4 e 2ie ", (A.24) 

rf n5 = 15 5 6 - 18i5V 9 " +21<5 4 e 2ie ", (A.25) 

d n6 = 35 e + 7i5 5 e l9n , (A.26) 

rf n7 = -<5 6 . (A.27) 

It is clear that the coefficients d mn and c mn have a very similar structure. The integral 
J i3 (/c) can then be rewritten as 

12 1.4 3 7 



18 H 2 M Pl p 3 

— ikrie+ikTn 



k T n 



' j fc/fco + « fc r n 

which can be integrated, as earlier, to yield 

12 1.4 3 7 



(A.28) 



x 



i k 



m=0 



r ( m, — — h ikr n ] —V (to, —ikr] c + ik r n ) 



(A.29) 
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The rest of the calculation in arriving at the final expression (160) proceeds just as in 
the case of J 13 . 

Appendix A. 2. The coefficients f n and g n 

The coefficients /„ are given by 

fo =1, (A.30) 

/i=-<+ 6 -f, (A.31) 

h fee - ' ( A - 66 > 

f - 3tp3 39 ^ (AM) 

h - "p (A.34) 

15p 6 27*p 9 

A = + % (A.37) 

h = (A.39) 

U = ^, (A.40) 
while g n are given by 

<? =3, (A.41) 
18 jo 3 

5i =-* + -p^> ( A - 42 ) 
12ip 3 27p 6 

^="^ + ^ (A.43) 

3 P 3 27i P 6 / a AA , 

3 *P 3 9 P 6 /A „ C X 

^ 4 = "# _ (A ' 45) 

95 = _ (A - 46) 

3 P 6 27 *P 9 

^ 6 = -IfcT + l/"' (A ' 47) 

3 ^P 6 9 P 9 /A .ON 

^ = —jjr - -yr, (A.48) 

= (A-49) 
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99 = 3 ~f, (A.50) 
9io = (A.51) 
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